Semicond. Sci. Technol. 19 (2004) R47–R58 PII: S0268-1242(04)75819-8
TOPICAL REVIEW
An efficient technique for solving the
arbitrarily multilayered electrostatic
problems with singularity arising from
a degenerate boundary
Shiang-Woei Chyuan
1, Yunn-Shiuan Liao
1and Jeng-Tzong Chen
2 1Department of Mechanical Engineering, National Taiwan University, Taipei, Taiwan, Republic of China2Department of Harbor and River Engineering, National Taiwan Ocean University, Keelung, Taiwan, Republic of China
E-mail: yeaing@iris.seed.net.tw Received 6 February 2004 Published 30 July 2004 Online atstacks.iop.org/SST/19/R47 doi:10.1088/0268-1242/19/9/R02 Abstract
Engineers usually adopt multilayered design for semiconductor and electron devices, and an accurate electrostatic analysis is indispensable in the design stage. For variable design of electron devices, the BEM has become a better method than the domain-type FEM because BEM can provide a complete solution in terms of boundary values only, with substantial saving in modelling effort. Since dual BEM still has some advantages over
conventional BEM for singularity arising from a degenerate boundary, the dual BEM accompanied by subregion technology, instead of tedious calculation of Fourier–Bessel transforms for the spatial Green’s functions, was used to efficiently simulate the electric effect of diverse ratios of permittivity between arbitrarily multilayered domain and the fringing field around the edge of conductors. Results show that different ratios of permittivity will affect the electric field seriously, and the values of surface charge density on the edge of conductors are much higher than those on the middle part because of fringing effect. In addition, if using the DBEM to model the fringing field around the edge of conductors, the minimum allowable data of dielectric strength for keeping off dielectric breakdown can be obtained very efficiently.
(Some figures in this article are in colour only in the electronic version)
1. Introduction
Electrostatics, as used here, is the term given to the study of the interactions between electrically charged bodies. Severely speaking electrostatics includes the study of the forces which hold individual electrons and ions together to form atoms, the chemical forces which bind the atoms together and much of the science of plasma physics and
electrochemistry [1]. When device dimensions are much less than the wavelength of electromagnetic radiation at a particular frequency, then the response of the system at that frequency can be considered quasistatic in that the emission, transmission, or absorption of electromagnetic radiation can be ignored [2]. For many microelectromechanical systems (MEMS) and electron devices, knowledge of the electric potential (V ) and
electric field intensity (E) are very important [3]. Therefore, electrical engineers are familiar with electrostatic problems, and diverse numerical methods have been regularly used in
MEMS and EM (electromagnetics) [4]. Among different
numerical approaches, finite element method (FEM), which is based on the representation and approximate solution of boundary value problems of engineering mathematics in terms of partial differential equations [5], and boundary element method (BEM) based on integral equations [6] have moved from being research tools for scientists to become powerful design tools for engineers. One of the main advantages of BEM, when compared to FEM, is that discretizations are restricted only to the boundaries, making data generation much
easier. The BEM is also ideally suited to the analysis of
external problems where domains extend to infinity, since discretizations are confined to the internal boundaries with no need to truncate the domain at a finite distance and impose artificial boundary conditions, and to problems involving some form of discontinuity or singularity, due to the use of singular fundamental solutions as test functions. Especially for variable design of electron devices, many laborious works of finite element modelling compared to those of boundary element model are needed because BEM can provide a complete solution in terms of boundary values only, with substantial saving in modelling effort. Therefore, there is no doubt that BEM has become a very appealing approach in numerical simulation of EM and MEMS [7–12] even if many engineers still use commercial package and waste much time to set up adverse FEM models in the design stage nowadays.
Because modern MEMS and electron devices design usually contains very thin conducting plates (e.g. a parallel-plate capacitor), the singularity problems arising from a degenerate boundary (The degenerate boundary refers to a boundary, two portions of which approach each other such that the exterior region between the two portions becomes infinitely thin.) are frequently formed, and it is well known that the coincidence of the boundaries gives rise to an ill-conditioned problem. The sub-domain technique in conventional BEM with artificial boundaries for degenerate boundary has been introduced to ensure a unique solution. The main drawback of the technique is that the deployment of artificial boundaries is arbitrary and, thus, cannot be implemented easily into an
automatic procedure. In addition, model creation is more
troublesome than in the single domain approach. To tackle such degenerate boundary electrostatic problems, dual BEM (DBEM) has been proposed in [13], and the above-mentioned boundary value problems can be solved efficiently in the original single domain if using DBEM.
The paper is organized as follows. Section 2involves
the reviews of electrostatic Green’s functions in layered dielectrics. In section3, we briefly introduce the procedure of dual boundary integral equation for electrostatic problems.
Numerical results are provided and compared in section 4
to establish the suitability and accuracy of the DBEM. Some remarks based on the reported results were discussed in section5. Finally, there is a concise conclusion in section6.
2. Reviews of electrostatic Green’s functions in layered dielectrics
Generally the multilayered medium structure has been recently developed and designed in many microwave circuits to provide
high-cost performance as commercial products, accurate and expedient prediction of electric performance of high-speed interconnects is critically dependent on efficient simulation of capacitance and inductance of layered dielectric substrates. While using the integral equation formulations formulated by the boundary integral equation techniques for solving electrostatic problems, how to search for a Green’s function is seemed to be the first and most important thing.
For the case of layered substrates, commonly utilized techniques use one of two types of Green’s functions. The first one is used as kernel of the integral equation (the free-space Green function), in which case, polarization charges on the dielectric interfaces need to be introduced as unknowns [14, 15]. In [14, 15], the spatial domain approach using the free-space Green function has been developed to calculate the capacitance and inductance matrices of multiconductor transmission lines located arbitrarily in a multilayered medium of finite extent. The potential in the medium is expressed in terms of the free charge at the conductor–dielectric interfaces and the total charge at the dielectric–dielectric interfaces. The second one is called the spatial Green function which is much more complicated than the free-space Green function because it is consistent with the boundary conditions at the dielectric interfaces [16–19]. Comparing the free-space Green function with spatial Green’s function, one can see that: (1) although the free-space Green’s function is simple, the number of unknowns involved in the discrete problem is large. (2) For spatial Green’s function, the only polarization charges required as additional unknowns are those at any finite dielectric boundaries and the number of unknowns in this case is reduced significantly. (3) Even though closed-form expressions in terms of an infinite series are available for the cases of single-layer microstrip and stripline configurations, the calculation of spatial Green’s functions for multilayered substrates requires the computationally tedious calculation of Fourier–Bessel transforms because it is necessary to relate the Fourier coefficients in each layer to the boundary conditions on the interfaces of the dielectric layers.
Because there are also many difficulties for irregular
geometrical design while using the afore-mentioned
conventional free-space and spatial Green’s functions, a systematic and robust methodology is needed for solving the arbitrarily multilayered electrostatic problems. However, many difficulties will still be faced if using DBEM for the multilayered domain because we also need more external boundary conditions at the interface between two dielectric media, which should be calculated. Since electromagnetic problems often involve media with different physical properties and require knowledge of the relation of the field quantities at the interface between two media, the DBEM accompanied by subregion technique was successfully used to determinate how the E and D (D= εE, electric flux density) vectors change in crossing the interface in this paper. For electrical engineering practices, since the numbers of elements and nodes for FEM are much higher than those of DBEM to get a reasonable result, the present DBEM accompanied by the subregion technique has great potential in multilayered electrostatic problems for future industrial applications.
3. Dual boundary integral equations
For a homogeneous medium, the governing equation of electrostatics can be written in the following form,
∇2V = −ρ/ε (1)
where∇2is the Laplacian operator. Equation (1) is known as
Poisson’s equation; it states that the divergence of the gradient of electric potential (V ) equals−ρ/ε for a simple medium, where ε is the permittivity of the medium and ρ is the volume density of free charges. At points in a simple medium where there is no free charge, equation (1) is reduced to
∇2V = 0 (2)
which is known as Laplace’s equation. Equation (2) plays a very important role in MEMS and EM. It is the governing equation for electrostatic problems involving a set of conductors, such as capacitors, maintained at different potentials. Once V is found from equation (2), E can be
determined from −∇V , and the charge distribution on the
conductor surfaces can be determined from surface charge density ρs= εEn.
The electrostatic problem consists of finding the unknown potential function (or V ) in the partial differential equation. In addition to the fact that satisfies equation (2) within a prescribed solution region , the potential function must satisfy certain conditions on B which is the boundary of . Usually these boundary conditions are the Dirichlet and Neumann types. Therefore, the governing equation and boundary conditions of electrostatic problems could be written in the following form.
Governing equation:
∇2(x)= 0, xin (3)
Dirichlet boundary condition:
(x)= f (x), xon B (4)
Neumann boundary condition:
∂(x)/∂nx = g(x), xon B, (5)
where f (x) and g(x) denote the known boundary data, and nx is the unit outer normal vector at the point x on the boundary. Based on the dual boundary integral equation formulation for electrostatic problem, we have
α(x)= CPV B T (s, x)(s)dB(s) − RPV B U (s, x)[∂(s)/∂ns] dB(s) (6) α[∂(x)/∂nx]= HPV B M(s, x)(s)dB(s) − CPV B L(s, x)[∂(s)/∂ns] dB(s), (7)
where the kernel functions, U (s, x) = ln(r), T (s, x) =
∂U (s, x)/∂ns, L(s, x) = ∂U(s, x)/∂nx, M(s, x) =
∂2U (s, x)/∂n
x∂ns, r= |s − x|, s and x being the position vectors of the points s and x, respectively, and nsis the unit outer normal vector at point s on the boundary (see figure1). In addition, RPV is the Riemann principal value, CPV is the Cauchy principal value, HPV is the Hadamard principal
Figure 1. Boundary element discretization for degenerate boundary
and non-degenerate boundary.
Figure 2. A figure sketch of the multi-domain.
Figure 3. An interface between two media.
value, and α depends on the collocation point (α = 2π for
an interior point, α = π for a smooth boundary, α = 0 for
an exterior point). The commutativity property of the trace operator and the normal derivative operator provides us with alternative ways to calculate the Hadamard principal value
analytically [12]. First, L’Hospital’s rule is employed in the limiting process. Second, the normal derivative of the Cauchy principal value should be taken carefully by using Leibnitz’ rule, and then the finite part can be obtained. The finite part has been termed the Hadamard principal
value or Mangler’s principal value. In the derivation of
dual equations, two alternatives can be applied to determine the Hadamard principal value as presented in appendix A. Generally, equation (6) is called singular boundary integral equation, and equation (7) is called hypersingular boundary integral equation. Since the hypersingular equation plays an important role in the degenerate problems, many researchers
have paid much attention to this. After discretizing the
boundary into 2N boundary elements, equations (6) and (7) are reduced to
[U ]2N×2N{t}2N×1= [T ]2N×2N{u}2N×1 (8)
[L]2N×2N{t}2N×1= [M]2N×2N{u}2N×1 (9)
where [U ], [T ], [L] and [M ] are the four influence matrices, {u} and {t} are the boundary data for the primary and the secondary boundary variables:{} and {∂/∂n}, respectively. These above-mentioned four influence matrices for interior and exterior fields are presented in appendix B. Generally for electrostatic problems without degenerate scale, the aforementioned influence matrix [U ] is nonsingular, either equation (8) or equation (9) can be solved by Gaussian elimination and LU decomposition very well. But for degenerate scale problem, [U ] matrix is singular and the rank is deficient, then the SVD (singular value decomposition) technique need to be used [20]. For multilayered electrostatic problems, the DBEM accompanied by subregion technique for multi-domain is needed, and the methods are presented in appendix C.
4. DBEM simulation of multilayered electrostatic problems
In order to demonstrate the suitability and efficiency of DBEM presented in this paper, two electrostatic problems were used. Basically for the nondegenerate boundary problems like the following case 1, either singular BEM ( just using the first kind kernels of U(s, x) and T(s, x)), hypersingular BEM (the second kind kernels of L(s, x), M(s, x)) or DBEM (using the first kind kernels of U(s, x), T(s, x) and the second kind kernels of L(s, x), M(s, x) in the meantime) accompanied by subregion technology can be used. But for the multilayered electrostatic problems with degenerate boundary like the following case 2, the DBEM accompanied by subregion technology plays an important role, and conventional BEM (singular BEM or hypersingular BEM) without external artificial boundaries cannot be used.
Case 1. Consider the potential problem shown in figure 4. The potentials at x= 0, x = a and y = 0 sides are zero while the potential at y = b side is V0. Try to find the potential
distribution in the multilayered domain.
Because it is not easy to obtain analytical data, FEM was used to calculate the reference data [21]. For convenience, the values of a, b and c are assumed to be 30, 20 and 10 µm. respectively. Four points will be analysed using rough mesh discretization (72 elements and 72 nodes, see
Figure 4. Figure of case 1 (without degenerate boundary).
Figure 5. The related DBEM mesh discretization of case 1 (without
degenerate boundary).
Figure 6. The related FEM mesh discretization of case 1.
figure 5) of conventional DBEM (or BEM) accompanied
by subregion technology, and compared with reference data computed from a larger FEM model (600 elements and 651 nodes, see figure 6). The results of electric potential under diverse numerical methods were listed in tables1and
2 and shown in figures7 and8. Comparing the results of
electric potential field (equipotential lines) using DBEM and FEM (see figures7and8), one can see that the difference of electric potential distribution is very little. Therefore, DBEM (or BEM) accompanied by subregion technology used in this
Figure 7. The results of electric potential (V ) (equipotential lines; ranging from +V0at top (red in online edition) to 0 at bottom (dark blue online)) of case 1 (R= 10) if using DBEM (left part) and FEM (right part).
Figure 8. The results of electric potential (V ) (equipotential lines; ranging from +V0at top (red in online edition) to 0 at bottom (dark blue
online)) of case 1 (R= 0.1) if using DBEM (left part) and FEM (right part).
Table 1. The results of electric potential (V ) of case 1 (R= 10) if using the conventional BEM, DBEM and FEM.
Results from Results from Results from Difference
Locations singular BEM hypersingular BEM DBEM Results from between
(x, y) U, T kernels L, M kernels U, T, L, M kernels FEM DBEM and FEM
(18.0, 3.00) 0.173 02V0 0.184 90V0 0.173 02V0 0.174 3103V0 −0.74% (4.00, 9.00) 0.274 48V0 0.283 12V0 0.274 48V0 0.280 9692V0 −2.31% (25.0, 16.0) 0.596 07V0 0.637 56V0 0.596 07V0 0.600 0305V0 −0.67% (5.00, 17.0) 0.674 92V0 0.712 73V0 0.674 92V0 0.679 071V0 −0.61%
Table 2. The results of electric potential (V ) of case 1 (R= 0.1) if using the conventional BEM, DBEM and FEM.
Results from Results from Results from Difference
Locations singular BEM hypersingular BEM DBEM Results from between
(x, y) U, T kernels L, M kernels U, T, L, M kernels FEM DBEM and FEM
(18.0, 3.00) 0.017 302V0 0.018 490V0 0.017 302V0 0.017 419 43V0 0.67% (4.00, 9.00) 0.027 448V0 0.028 312V0 0.027 448V0 0.028 100 59V0 −2.32% (25.0, 16.0) 0.480 640V0 0.511 030V0 0.480 640V0 0.488 333 13V0 −1.58% (5.00, 17.0) 0.589 690V0 0.615 380V0 0.589 690V0 0.592 9200V0 −0.54%
paper seems to be a very efficiently numerical method for the multilayered electrostatic problems without degenerate boundary (e.g. thin strip conductors). From the results of electric potential V(x, y) listed in tables 1 and2, one can see that the differences between FEM and DBEM are lower than 2.4%. In addition, both singular BEM and hypersingular BEM accompanied by subregion technology can deal with this nondegenerate boundary problem well.
In order to investigate the effect of diverse ratios of permittivity (R = ε1/ε2)between subdomains 1 and 2, the
results of electric potentials of Y= 10 µm interface line and
X= 15 µm lines from DBEM are shown in figures 9 and
10. From figures9and10, one can see that the values of R seriously affect the distribution of voltages and engineers can adjust R according to their needs.
Case 2. Similar to case 1, but insert two thin strip conductors
(length = 10 µm) into the centre part of these two media
shown in figure4. If the voltage of the upper strip conductor is 0 and that of the lower one is V0, determine the electric field
and charge distributions.
Because these two strip conductors form a new degenerate boundary and a closed contour formed by two lines C+and C−
(see figure1) for the degenerate boundary, the singular BEM or hypersingular BEM without external artificial interface cannot cope with this degenerate case, but the electric potential can be calculated just in one run if using DBEM. Similar to case 1, FEM was also used to calculate the reference data, but the linear shape function of finite element model shown in figure6was replaced by parabolic shape function in order to model the fringing effect near the edge of these two strip
0 0.15 0.3 0.45 0.6 0.75 0 5 10 15 20 25 30 X Coordinate ( m) E lectr ic Po tentia ls (V ) R=10. R=5. R=2. R=1. R=0.5 R=0.2 R=0.1 µ
Figure 9. Results of electric potentials (V ) of Y= 10 µm interface line under diverse ratios of permittivity (R) between subdomains 1 and
2—case 1 (unit: V0). 0 0.175 0.35 0.525 0.7 0.875 1.05 0 4 8 12 16 20 Y Coordinate ( m) E lectric P o tentia ls ( V ) R=10. R=5. R=2. R=1. R=0.5 R=0.2 R=0.1 µ
Figure 10. Results of electric potentials (V ) of X= 15 µm line under diverse ratios of permittivity (R) between subdomains 1 and 2—case 1
(unit: V0).
conductors. Four points will also be analysed using more dense mesh discretization (180 elements and 156 nodes, see
figure11) of DBEM accompanied by subregion technology,
and compared with reference data computed from a larger FEM model (600 elements and 1901 nodes). The results of electric potential under diverse numerical methods were listed
in table 3 and shown in figure 12. Comparing the results
of electric potential field of DBEM with those of FEM (see figure12), one can see that the difference of electric potential distribution is also very little and the fringing effect near the edge of these two strip conductors can be modelled very well. From the results of electric potential V(x, y) listed in table3 and shown in figure12, one can find that (1) the conventional BEM without external artificial boundaries cannot solve the electrostatic problem with degenerate boundaries such as this case. (2) DBEM accompanied by subregion technology can analyse the multilayered electrostatic problem with degenerate boundary very efficiently. (3) The fringing effect near the edge of these two strip conductors can be simulated very well using DBEM. (4) The differences between FEM and DBEM are lower than 5.0%. In addition, the results of electric potentials
of Y= 10 µm interface line and X = 15 µm lines from DBEM
are also shown in figures13and14in order to investigate the effect of diverse ratios of permittivity (R) between subdomains 1 and 2. From figures13and14, one can see that the values of R also critically influence on the distribution of electric
Figure 11. The related DBEM mesh discretization of case 2 (with
degenerate boundary).
potentials like case 1 (see figures9and10). Comparing with the results of nondegenerate case, the electric potential field of degenerate case is more complicated because of degenerate boundary.
Besides electric potentials, the distributions of normal electric field (En) on the top and bottom sides of upper
conductor under diverse R can be shown in figures 15 and
Figure 12. The results of electric potential (V ) (equipotential lines; ranging from +V0at top and lower middle area (red in online edition) to 0 at bottom and upper middle area (dark blue online)) of case 2 (R= 10) if using DBEM (left part) and FEM (right part).
0 0.15 0.3 0.45 0.6 0.75 0.9 0 5 10 15 20 25 30 X Coordinate ( m) E lectr ic Po ten tia ls (V ) R=10. R=5. R=2. R=1. R=0.5 R=0.2 R=0.1 µ
Figure 13. Results of electric potentials (V ) of Y= 10 µm interface line under diverse ratios of permittivity (R) between subdomains 1 and
2—case 2 (unit: V0). 0 0.175 0.35 0.525 0.7 0.875 1.05 0 4 8 12 16 20 Y Coordinate ( µm) E lectric P o tentia ls ( V ) R=10. R=5. R=2. R=1. R=0.5 R=0.2 R=0.1
Figure 14. Results of electric potentials (V ) of X= 15 µm line under diverse ratios of permittivity (R) between subdomains 1 and 2—case 2
(unit: V0).
Table 3. The results of electric potential (V ) of case 2 (R= 10) if using the conventional BEM, DBEM and FEM.
Results from Results from Results from Difference
Locations singular BEM hypersingular BEM DBEM Results from between
(x, y) U, T kernels L, M kernels U, T, L, M kernels FEM DBEM and FEM
(24.0, 16.5) NA NA 0.521 81V0 0.514 489 8V0 1.42%
(6.50, 12.0) NA NA 0.238 01V0 0.230 1575V0 3.41%
(22.5, 6.00) NA NA 0.346 38V0 0.363 8855V0 −4.81%
(4.00, 3.50) NA NA 0.106 23V0 0.110 8643V0 −4.18%
bottom side of upper conductor are obviously dependent on the values of R and the location to the left corner of conductor (locx), and the smaller the R is, the larger the En,2is. Unlike
En,2, the values of En,1 on the top side are only apparently counting on the value of locx, and the effect of the values of
R can be ignored (see figure 15). Similar to the results of upper conductor, the distributions of Enon the top and bottom sides of lower conductor under diverse R can also be shown in figures 17and18. From figure17, one can see that the values of En,3on the top side of lower conductor are obviously
0 0.1 0.2 0.3 0.4 0.5 0 2 4 6 8 10
Location to the left side of conductor ( m)
No rm a l electr ic field (E n ) R=0.1 R=0.2 R=0.5 R=1. R=2. R=5. R=10. µ
Figure 15. The distribution of normal electric field (En) on the top side of upper conductor under diverse ratios of permittivity (R) between
subdomains 1 and 2 (unit: V0/µm).
0.00 0.10 0.20 0.30 0.40 0 2 4 6 8 10
Location to the left side of conductor ( m)
No rm a l electr ic field (E n ) R=0.1 R=0.2 R=0.5 R=1. R=2. R=5. R=10. µ
Figure 16. The distribution of normal electric field (En) on the bottom side of upper conductor under diverse ratios of permittivity (R)
between subdomains 1 and 2 (unit: V0/µm).
-0.55 -0.44 -0.33 -0.22 -0.11 0 0 2 4 6 8 10
Location to the left side of conductor ( m)
No rm a l electr ic field (E n ) R=0.1 R=0.2 R=0.5 R=1. R=2. R=5. R=10. µ
Figure 17. The distribution of normal electric field (En) on the top side of lower conductor under diverse ratios of permittivity (R) between
subdomains 1 and 2 (unit: V0/µm).
-0.55 -0.44 -0.33 -0.22 -0.11 0 0 2 4 6 8 10
Location to the left side of conductor ( m)
N o rm a l el ect ri c fi el d (E n) R=0.1 R=0.2 R=0.5 R=1. R=2. R=5. R=10. µ
Figure 18. The distribution of normal electric field (En) on the bottom side of lower conductor under diverse ratios of permittivity (R)
between subdomains 1 and 2 (unit: V0/µm).
dependent on the values of R and the locx, but the values of En,4 on the bottom side are only apparently counting on the value of
locx, and the effect of the values of R can be also neglected (see figure18). From figures15to18, one can find that different
-1 -0.8 -0.6 -0.4 -0.2 0 0 2 4 6 8 10
Location to the left side of conductor ( m)
S u rfac e c h ar ge de n si ty R=0.1 R=0.2 R=0.5 R=1. R=2. R=5. R=10. µ
Figure 19. The distribution of surface charge density (ρs) of upper conductor under diverse ratios of permittivity (R) between subdomains
1 and 2 (unit: ε1V0/µm2). 0 0.22 0.44 0.66 0.88 1.1 0 2 4 6 8 10
Location to the left side of conductor ( m)
S u rfac e c h ar ge de n si ty R=0.1 R=0.2 R=0.5 R=1. R=2. R=5. R=10. µ
Figure 20. The distribution of surface charge density (ρs) of lower conductor under diverse ratios of permottivity (R) between subdomains
1 and 2 (unit: ε2V0/µm2). -6 -4 -2 0 2 4 6 0 2 4 6 8 10
Ratios of permittivity between Subdomain 1 and Subdomain 2 (R) C h ar ge (Q ) Upper conductor Lower conductor
Figure 21. The distribution of charge (Q) on both upper and lower conductors under diverse ratios of permittivity (R) between subdomains
1 and 2 (unit: upper conductor: ε1wV0; lower conductor: ε2wV0).
R will affect the electric field seriously, and the values of En on the edge of conductors are much higher than those on the middle part because of the fringing effect.
As we know the value of R can play a very important role, let us investigate the effect of R for charge distribution. Because the charge distribution on the conductor surfaces can be determined from ρs = εEn(the normal component of the electric field Enat a conductor boundary is equal to the surface charge density ρson the conductor divided by the permittivity
ε [1, 2]) if ε is a constant. From figures19and20, we can see that the values of surface charge density (ρs) on the both upper and lower conductors are obviously dependent on the value of locx, but not very obvious for the value of R. As we all know that the fringing effect around the edge of inserted conductors is as clear as day like En. From these results, we can see the ρsat the edges becomes much larger than that at
the centre for conductors because of fringing effect. If the width of conductor is w, the distribution of charge (Q) on both upper and lower conductors under diverse R can be shown in figure21. From figure21, we can see that the larger the value of R is, the larger the Q of lower capacitor is, but the larger the value of R is, the less the Q of upper capacitor is.
5. Discussions
(1) For the electrostatic problems without degenerate boundary like case 1, singular BEM, hypersingular BEM and DBEM are all very usefully numerical tools. For the electrostatic problems with singularity arising from degenerate boundary such as case 2, it is well known that the coincidence of the boundaries will give rise to
an ill-conditioned problem if using conventional BEM. Therefore DBEM presented in this paper became a very efficient solver for modelling the fringing field around the edge of electron devices and the conventional BEM without external artificial boundaries cannot be used for the electrostatic problem with degenerate boundaries. (2) By way of DBEM, one can see that the values of normal
electric field (En) on the bottom side of upper and top side of lower conductors are obviously dependent on the values of ratios of permittivity (R) between subdomains 1 and 2 and the location to the left corner of conductor (locx), but the values of Enon the top side of upper and bottom sides of lower conductors are only apparently counting on the value of locx, and the effect of the values of R can be neglected. Results also show that different R will affect the electric potential (V ) and field seriously, and the values of electric field and surface charge density (ρs) on the edge of conductors are much higher than those on the middle part because of fringing effect.
(3) As we know from the theory of electrostatics [1, 2], the electrons will be pulled out of the molecules completely if the E is very strong, and the electrons will accelerate under the influence of the E, collide violently with the molecular lattice structure, then cause permanent dislocations and damage in the material. While avalanche effect of ionization due to collisions occurs, the material will become conducting, and large currents (I ) may result. This phenomenon is called a dielectric breakdown, and the maximum electric field intensity that a dielectric material can withstand without breakdown is the dielectric strength of the material (e.g. for air is 3 × 106 V m−1). So
the accurate and efficient modelling for fringing effect around the edge of electron devices is very important for performance because we need to know the minimum allowable data of dielectric strength for keeping off
dielectric breakdown. Therefore, we recommend our
DBEM here because it is more simple and efficient than other conventional BEMs.
(4) For the arbitrarily multilayered electrostatic problems with singularity arising from degenerate boundary, DBEM accompanied by subregion technology can efficiently model the physical behaviour of interface between subdomains 1 and 2. In addition, for variable design of electron devices, the DBEM has become a better method than the domain-type FEM because DBEM can provide a complete solution in terms of boundary values only, with substantial saving in modelling effort.
Especially for calculating the ρs and charge (Q) of
conductors according to diverse design, only the part of Enon the surface of conductors needed to be known, then the DBEM can work very well.
6. Conclusions
The dual boundary integral formulation accompanied by subregion technology for arbitrarily multilayered electrostatic problems has been presented in this article. Comparisons of the electric potentials between the FEM and DBEM analyses were discussed with respect to degenerate and non-degenerate designs in order to demonstrate the suitability and efficiency
of computational method presented in this paper. It has been shown that the DBEM accompanied by subregion technology in the context of the present formulation is particularly suitable for the multilayered electrostatic problem with singularity arising from a degenerate boundary. For electrical engineering practices, since the major effort is model creation, the present DBEM, free from the development of an artificial boundary, has great potential for industrial applications.
Appendix A
In the derivation of dual equations, two alternatives can be applied to determine the Hadamard principal value as follows [12]:
(1) Trace operator first and differential operator second HPV B M(s, x)u(s)dB(s) = ∂ ∂nx CPV B T (s, x)u(s)dB(s) . (A1)
For simplicity, constant element is adopted, i.e. u(s)= 1, and (A1) reduces to d dx CPV c a −1 (x− s)ds = d dx x−ε a −1 (x− s)ds + c x+ε −1 (x− s)ds = 1 a− x − 1 c− x after using the Leibnitz rule.
(2) The differential operator first and trace operator second HPV B M(s, x)u(s)dB(s)= lim y→x B M(s, y)u(s)dB(s). (A2) Similarly, constant element scheme can simplify (A2) into
lim y→0 c a 1 (x− s)2+ y2ds = lim y→0 1 y tan−1 y a− x − tan−1 y c− x = 1 a− x − 1 c− x after using lim α→0 tan−1α α = 1 and tan−1α+ tan−11 α = π 2. Appendix B
In dual BEM, the linear algebraic equations for an interior electrostatic problem discretized from the dual boundary integral equations can be written as [12, 13]
Tpqi {q} = Upqi {∂/∂n}q (B1) Mi pq {q} = Li pq {∂/∂n}q (B2)
where {q} and {∂/∂n}q are the boundary potential and flux, and the subscripts p and q correspond to the labels of the collocation point and integration element, respectively. For the exterior electrostatic problem, we have [12, 13]
Tpqt {q} = Upqe {∂/∂n}q (B3) Mpqe {q} = Lepq {∂/∂n}q. (B4)
The influence coefficients of the four square matrices [U ], [T ], [L] and [M ] can be represented as
Upq = RPV Bq U (sq, xp)dB(sq) (B5) Tpq = −πδpq+ CPV Bq T (sq, xp)dB(sq) (B6) Lpq = πδpq+ CPV Bq L(sq, xp)dB(sq) (B7) Mpq= HPV Bq M(sq, xp)dB(sq) (B8)
where Bq denotes the qth element and δpq = 1 if p = q,
otherwise it is zero. The explicit form will be derived in the following section. According to the dependence of the out-normal vectors in these four kernel functions for the interior and exterior electrostatic problems, their relationship can be easily found: Upqi = Upqe (B9) Mpqi = Mpqe (B10) Tpqi = −Tpqe if p= q; Tpqi = Tpqe if p= q. (B11) Lipq= −Lepq if p= q; Lipq= Lepq if p= q. (B12) Appendix C
(1) Subregion technology of BEM. Before studying the physical behaviour of interface between two media of electrostatic problems, the concept of subregion technology for singular BEM ([U ] and [T ]) or hypersingular BEM ([L] and [M ]) was introduced first. From figure2, a figure sketch of the multi-domain could be found, whereu1
c
is the boundary data of degenerate boundary (C+and C−), andu1
f ,u2 f ,t1 f ,
tf2 are the unknowns of the interface between subdomains 1 and 2. Since the degenerate boundary on C+and C−as shown
in figure 2results in double unknowns, Equation (8) or (9) can provide an additional equation by collocating the point x on C+and C−. By dividing the domain into two subdomains
(subdomains 1 and 2) and using singular BEM or hypersingular BEM influence matrices for each subdomain, we have the two equations from equation (8) or (9) as follows:
T1 cc T 1 cf T1 f c Tff1 u1 c u1f = U1 cc U 1 cf U1 f c Uff1 t1 c tf1 (C1) and Tcc2 Tcf2 T2 f c Tff2 u2 c u2f = Ucc2 Ucf2 U2 f c Uff2 t2 c tf2 (C2)
where the superscripts 1 and 2 are the labels of the subdomains and the subscripts c and f denote the complementary and interface sets for u and t, respectively. Since the unknown pairs ofu1f ,u2f ,tf1 andtf2 are introduced in the artificial boundary as shown in figure2, two constraints of the continuity and equilibrium conditionsu1
f =u2 f andt1 f = −t2 f are necessary for the interface if ε1= ε2. By assembling (C1)
and (C2) with two constraints of the continuity and equilibrium conditions for the interface, we have
Ucc1 Ucf1 0 U1 f c Uff1 0 0 −U2 cf Ucc2 0 −U2 ff Uf c2 t1 c t1 f tc2 = Tcc1 Tcf1 0 T1 f c Tff1 0 0 T2 cf Tcc2 0 T2 ff Tf c2 u1 c u1 f u2c . (C3) After solving (C3), all the unknownst1
c ,t2 c ,t1 f andu1 f can be obtained, and u2
f
andt2
f
also can be calculated from the two constraints of the continuity and equilibrium conditions for the interface.
(2) Boundary conditions for electrostatic fields of the interface between two media. If there is an interface between two
general media shown in figure 3, the normal component of
D field is discontinuous across an interface where a surface charge exists—the amount of discontinuity being equal to the surface charge density. When two dielectrics are in contact with no free charges at the interface, ρs= 0, we have
ε1E1n= ε2E2n. (C4)
In addition, the electric potential and the tangential component of an E field are all continuous across an interface [1, 2]. (3) Subregion technology for multilayered electrostatic problems. For multilayered domain (ε1 = ε2) shown in
figure3, by assembling (C1) and (C2) with two constraints of the continuity and boundary conditions (C4) for the interface
u1f =u2f andtf1 = −Rtf2 , where R = ε1/ε2is the
ratio of permittivity between subdomains 1 and 2, we have U1 cc Ucf1 0 U1 f c Uff1 0 0 −RU2 cf Ucc2 0 −RUff2 Uf c2 t1 c t1 f tc2 = T1 cc Tcf1 0 T1 f c Tff1 0 0 T2 cf Tcc2 0 Tff2 Tf c2 u1 c u1 f u2c . (C5) After solving (C5), all the unknownst1
c ,t2 c ,tf1 andu1f can be obtained, andu2
f
andt2
f
also can be calculated from the two constraints of the continuity and boundary conditions for the interface.
References
[1] Cross J A 1987 Electrostatics: Principles, Problems and
Applications (Bristol: Adam Hilger)
[2] Haus H A and Melcher J R 1989 Elctromagnetic Fields and
Energy (Englewood Cliffs, NJ: Prentice Hall)
[3] Liu D C 2003 On the theory of the modal expansion method for the EM fields in side a closed volume J. Phys. D: Appl.
Phys. 36 1629–33
[4] Sadiku M N O 1992 Numerical Techniques in
Electromagnetics (Boca Raton, FL: CRC Press)
[5] Jin J 2002 The Finite Element Method in Electromagnetics (New York: Wiely)
[6] Aliabadi M H 2002 The Boundary Element Method (New York: Wiely)
[7] Senturia S D, Aluru N and White J 1997 Simulating the behavior of MEMS devices: computational methods and needs IEEE Comput. Sci. Eng. Mag. 30–43
[8] Wiersig J 2003 Boundary element method for resonances in dielectric microcavities J. Opt. A: Pure Appl. Opt. 5 53–60 [9] Ye W, Mukherjee S and MacDonald N C 1998 Optimal shape
design of an electrostatic comb drive in microelectromechanical systems IEEE J.
Microelectromech. Syst. 7 16–26
[10] Chakravorti S and Steinbigler H 2000 Boundary element studies on insulator shape and electric field around HV insulators with or without pollution IEEE Trans. Dielectr.
Electr. Insul. 7 167–76
[11] Chyuan S W, Liao Y S and Chen J T 2004 Computational study of variations in gap size for the electrostatic levitating force of MEMS device using dual BEM Microelectron. J. at press
[12] Liao Y S, Chyuan S W and Chen J T 2004 An alternatively efficient method (DBEM) for simulating the electrostatic field and levitating force of a MEMS combdrive
J. Micromech. Microeng. 14 1258–69
[13] Chyuan S W, Liao Y S and Chen J T 2003 An efficient method for solving electrostatic problems IEEE Comput. Sci. Eng.
Mag. 52–8
[14] Rao S M, Sarkar T K and Harrington R F 1984 The
electrostatic field of conducting bodies in multiple dielectric media IEEE Trans. Microwave Theory Tech. 32 1441–8 [15] Wei C, Harrington R F, Mautz J R and Sarkar T K 1984
Multiconductor transmission lines in multilayered dielectric media IEEE Trans. Microwave Theory Tech. 32 439–50 [16] Li K, Atsuki K and Hasegawa T 1997 General analytical
solutions of static Green’s functions for shielded and open arbitrarily multilayered media IEEE Trans. Microwave
Theory Tech. 45 2–8
[17] Cangellaris A and Yang L 2001 Rapid calculation of electrostatic Green’s functions in layered dielectrics IEEE
Trans. Magn. 37 3133–6
[18] Smith P M 2001 Dyadic Green’s functions for multi-layer SAW substrates IEEE Trans. Ultrason., Ferroelectr. Freq.
Control 48 171–9
[19] Tan E L 2002 Electrostatic Green’s functions for planar multilayered anisotropic media IEE Proc., Microw.
Antennas Propag. 149 78–83
[20] Chen J T and Lin S R 2002 On the rank-deficiency problems in boundary integral formulation using the Fredholm alternative theory and singular value decomposition technique Proc. 5th World Congress on Computational
Mechanics (Vienna, 7–12 July)
