Galerkin method for discs capacitors

Available online at www.sciencedirect.com

ScienceDirect Mathematics and Computers in Simulation 166 (2019) 365–381 www.elsevier.com/locate/matcom

Galerkin method for discs capacitors Giampiero Paffuti University of Pisa, Physics Department, Largo Pontecorvo 7, Pisa, Italy Received 1 February 2018; received in revised form 25 January 2019; accepted 16 June 2019 Available online 26 June 2019

Abstract The physics of nano and micro electro-mechanical systems would benefit from an accurate knowledge of electrostatics forces between electrodes. Quite complicated behaviors can arise and a precise computation of capacity matrices, from which forces can be derived, is required. We develop, for circular thick electrodes, a method based on the Galerkin approximation. Almost all the matrix elements needed are evaluated analytically and consequently high precision results can be obtained. In parallel we provide an implementation of the Boundary Elements Method (BEM) which gives a check of the computations and can be used independently. We test our method on the well studied problem of two flat discs capacitor, obtaining results in excellent agreement with expected asymptotic estimates. c 2019 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights ⃝ reserved. Keywords: Galerkin method; Coefficients of capacitance; Boundary Elements Method

1. Introduction In this work we propose a quite accurate numerical method for the calculation of the capacitance matrix for cylindrical capacitors, the method is also easily adaptable to the computation of the capacity of cylindrical shaped single conductors. We use a Galerkin procedure based on Gegenbauer and Jacobi polynomials. To check the procedure we give a Boundary Elements Method (BEM) integrated both on charge and field patches. The need for accuracy in this field stems from the following physical problem. It has been shown recently [20,22] that a repulsive electrostatic force is present for two close flat and equal electrodes, this has been verified both for discs [13] and squares [14]. It is natural to extend the study to more realistic tick electrodes, the method presented in this work is suitable for the case of circular conductors. In this work we do not discuss in detail the physical applications [15,21] of the method limiting ourselves to the introduction of general mathematical and numerical aspects of the problem, however we will present a numerical test on the system composed by two flat discs. The problem of the two (flat) discs capacitor has a long tradition in the physical and mathematical literature. In particular the behavior of capacitance for close approach of the plates has been studied by Kirchhoff in the pioneering work [8] and the results have been extended in [1,3,24,28]. The formulation of the problem in terms of a dual integral equation in [2,11,18,26] allowed a powerful analytical and numerical framework for the problem E-mail address: [email protected]. https://doi.org/10.1016/j.matcom.2019.06.009 c 2019 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights 0378-4754/⃝ reserved.

366

G. Paffuti / Mathematics and Computers in Simulation 166 (2019) 365–381

and permitted a rigorous derivation of Kirchhoff’s result [6]. A generalization of these results has been presented in [20,22] where the numerical accuracy has been improved and now agrees in the short distance region with the above asymptotic estimates. We will use this particular case as a test of our method in Section 4. For the problem of two thick conductors the only known result, to our best knowledge, is the Kirchhoff [8] calculation, which gives a part of the capacitance matrix at short distances, not sufficient for the computation of forces in the general case. The main tool in this work is the Galerkin method. The method is well known [25] and has been applied to systems similar to those studied in this work, see for example [7] for the case of a single flat disc with a dielectric substrate. The novelty in this work is the presentation of an accurate method for studying the whole capacitance matrix for a capacitor, built by two thick conductors. The accuracy achieved allows to study different physical properties, such as the forces, the quadrupole moments, the polarizability of these systems. To our best knowledge the only comparable framework is the case of a capacitor made of two flat disks, where the integral Love [11] equation can be used [20,22]. A crucial improvement of the present approach is to allow a thickness, making these models more realistic. New features appear in this new context as the qualitative change in the nature of the forces depending on the thickness. In this paper we focus on the mathematical aspects of the problem referring to [15,21] for a detailed discussion of the physical aspects. The paper is organized as follows. In Section 2 we describe the systems under study and define the relevant physical variables. In Section 3 we describe the application of the Galerkin method to our problem and present an implementation of BEM. Section 4 contains a numerical test on a two flat discs condenser while in Section 5 we collect the explicit expression of the matrices used in the Galerkin method. Section 6 contains some final remarks. 2. Description of the system In a system of conductors the charges Q i and the potentials Vi are related by the symmetric capacitance matrix Ci j Qi =

∑ j

Ci j V j ;

Vi =

∑

Mi j Q j

(1)

j

Ci j is the (symmetric) capacitance matrix, Mi j are the matrix element of the inverse. In terms of these matrices the energy takes the form 1∑ 1∑ Ci j Vi V j = Mi j Q i Q j (2) W = 2 ij 2 ij In this work we consider the case of two equal thick disks, i.e. two equal cylinder, at close contact with one base. The radius of the basis will be denoted by a and the thickness by b. The symmetry of the problem imply C11 = C22 and an analogous relation holds for Mi j . The distance between the nearest surfaces of the two electrodes will be denoted by ℓ. We will use the dimensionless variables κ = ℓ/a, τ ≡ 2h/a = b/a, i.e. a is our unit of length. The geometry is shown in Fig. 1, where the notations for a single cylinder are also given. In [13,14,20,22] it has been pointed out that the capacitance coefficients can be organized in a hierarchy according their behavior for κ → 0. In particular the sum C11 + C12 stays finite in the limit κ → 0 then it is useful to consider as independent quantities C11 − C12 ; C g1 = C11 + C12 (3) 2 C is the usual mutual capacity, i.e. the quotient between the charge on one electrode and the potential difference when the two electrodes are held at opposite charges. In terms of the above quantities the force between the two electrodes can be written as C=

F =−

∂ 1 (Q 1 + Q 2 )2 ∂ 1 (Q 1 − Q 2 )2 ∂ W = C g1 + C. 2 ∂ℓ 4 C g1 ∂ℓ 8 C2 ∂ℓ

(4)

In the works mentioned above it is shown that the second term in (4) produces a constant attractive force at short distances, for planar contacts. The first term is in general repulsive, and in the particular case of two equal (flat)

G. Paffuti / Mathematics and Computers in Simulation 166 (2019) 365–381

367

Fig. 1. Sketch of the systems studied in this work: a capacitor composed by two tick discs and a single cylinder. In Section 4 we consider as a numerical test the particular case of two flat discs, i.e. with b = 0.

discs it is logarithmically divergent [13,20] for κ → 0. These results follow from the known results for flat discs, valid for κ → 0: { [( ]} [ ( ) ]} { 1 κ )2 1 1 1 C →a κ log −2 , (5) + log 16π −1 +a 4κ 4π κ 16π 2 16π and [ C g1 → a

] 1 κ ( κ ) + 1 − log( ) . π 2π 2 π

(6)

For the discussion of these results and the comparison with numerical computations we refer to the above references, a summary is given in [20]. For thick discs, to our best knowledge, the only known result is the original computation of Kirchhoff [8,23], confirmed in [24]: { } ( ( 1 1 [ κ τ) τ) τ τ] C →a − 1 + log − 1+ log 1 + + log (7) 4κ 4π 16π κ κ κ κ and nothing is analytically known for C g1 . The methods of this work allow the computation of C and C g1 , in particular in [15,21] it is shown that Eq. (7) requires a finite correction depending on τ . The system we are considering has a plane of symmetry, orthogonal to the cylinder axis and in this case the study of the coefficients C, C g1 brings also some simplifications from the mathematical point of view. From (1) it follows that for the determination of the two coefficients it suffices to consider the case in which the conductors are at the same potential or at the opposite potential, having in the two cases Q 1 = (C11 + C12 )V = C g1 V ;

Q 1 = (C11 − C12 )V = 2C V

(8)

Our goal is then reduced to the computation of Q 1 . For a system with a symmetry plane the charge distribution at equilibrium is symmetric or antisymmetric with respect to this plane in the two cases, i.e. if on conductor 1 there is a charge density σ (x), on the reflected point x R on conductor 2 there is a charge density σ (x R ) = ±σ (x). The equilibrium charge in the two cases is then computed by solving two similar integral equations [ ] ∫ 1 1 ± σ (y) = V (9) |x − y R | y∈S1 |x − y| The variable y runs on the conductor 1, together with the reflected coordinate y R on the second conductor. Once computed σ on the conductor 1 its integration gives Q 1 and thus the two coefficients C g1 and C using (8). Let us consider the case under study, two equal discs, and choose for the upper cylinder the coordinate axis z centered in the body and positive upwards, with range −h ≤ z ≤ +h. In our configuration the density depends only on z and on the distance r from the symmetry axis. If we label the points of the second discs with z pointing

368

G. Paffuti / Mathematics and Computers in Simulation 166 (2019) 365–381

downward, in the same interval, and with the origin in medium plane of the second disc we have, by symmetry, for the two charge densities: σ (1) (r, z) = σ (r, z) ;

σ (2) (r, z) = ±σ (r, z)

i.e. the same function describes both densities. The explicit implementation of this procedure for the discs is given in Section 3. A similar and even simpler approach can be used for a single conductor, a single cylinder in our case. For a single conductor the determination of the superficial charge density allows not only the determination of the capacity, but also the determination of other electrostatic parameters, the most important of which is the (longitudinal) quadrupole moment for unit charge: ∫ 1 σ (2z 2 − (x 2 + y 2 )) (10) dzz = Q surf A slight modification of the procedure allows the computation of the longitudinal polarizability. The presence of a uniform electric field E directed along the z axis do not destroy the symmetry of the problem: the equilibrium charge density induced on the cylinder will be antisymmetric with respect to the median plane of the body and such to cancel the external potential. For a cylinder this means that the charge density σ must determine a constant potential ±Eh on the basis and E z on the lateral surface. Once found the charge density we can immediately found the longitudinal polarization by computing the induced dipole moment ∫ αzz E = zσ (11) surf

In the following we put for brevity α = αzz , D = dzz . Apart for their intrinsic importance these parameters are relevant for the large distance behavior of the capacitance coefficients [12], then their precise values can help to test the consistency of the procedure. For a single conductor we will use often the symbol L = 2h ≡ ax to denote the length of the cylinder, while for two conductors the reader has to be careful to distinguish between the distance between nearest surfaces, denoted by ℓ = κa as stated above, and the distance between the centers of the cylinders, denoted by d in the formulas below and in Section 5. Clearly d = ℓ + b ≡ ℓ + 2h . The explicit formulas for the matrix elements will be given in terms of d and h. 3. Numerical methods In this work we present an approximate solution of Eq. (9) for the capacitor composed by two thick discs with two different methods: the Galerkin method and a particular form of the Boundary Element Method (BEM). 3.1. Galerkin method Consider the electrostatic problem (9) and let us expand the unknown charge density in a complete basis of functions: ∞ ∑ σ (x) = cn f n (x) n=1

Inserting in (9) and performing the integral we transform the integral equation in an infinite set of equations ∫ ∞ ∑ An (x)cn = V ; An (x) = K (x, y) f n (y) (12) y

n=1

K stays for the integral (Coulombic) kernel and x is a generic point on the surface of the conductors. Projecting again the system (12) on the basis functions we convert the initial problem in an infinite set of linear equations ∫ ∫ ∑ Amn cn = Vm ; Amn = f m (x)An (x) Vm = f m (x)V (13) n

x

x

G. Paffuti / Mathematics and Computers in Simulation 166 (2019) 365–381

369

Except for very particular cases the system (13) does not admit an analytical known solution. The Galerkin approximation consists in a truncation of the system, i.e. we consider only functions with n ≤ N , in this way the problem reduces, in principle, to the solution of a system of N linear equations. In our case we have to determine three different functions: the surface densities on the two basis, σ B(1) (r ), σ B(2) (r ) which depend only on r , and the density on the lateral surface, σ L (z), which depends only on z. For the problem at hand we choose a slight generalization of the basis suggested in [4,27] ∑ ∑ σ L = aV an gn (z) ; σ B = aV bn f n (r ) ; (14) n

where ( 1( r 2 ) p−1 Γ (n) r2 ) (0, p−1) 1 − P 1 − 2 n−1 a2 a2 2 p−1 Γ (n − 1 + p) a2 ( ( ) s 2 )s−1/2 2 Γ (s)Γ (n) s z 1 z C gn (z) = 1− 2 2πah h Γ (n − 1 + 2s) n−1 h f n (r ) =

(15) (16)

(α,β)

Pn and Cns are respectively Jacobi and Gegenbauer polynomials. The coefficients for the two basis will be denoted (1) by bn , bn(2) while an are the coefficients for the lateral surface. The parameters s, p are in principle arbitrary, but clearly it is better to choose them in such a way the functions embody the expected singular behavior at the edges. For a solid cylinder the surfaces intersect at ninety degrees and the expected singularity is of the form ξ −1/3 , where ξ is the distance from the edges, see [9] (problem 3 of § 3) or [17]. This is realized by the choice s = 1/6, p = 2/3. For hollow cylinders only the Gegenbauer polynomials are needed while for discs only Jacobi polynomials are used. The edges of these systems require a charge density with a behavior ξ −1/2 , where ξ is the distance from the edge, see references [9,17] cited above, then for the two systems the natural choice for the parameters of the polynomials is s = 0, p = 1/2 respectively. In [21] the freedom in the choice of the parameters s, p studying the asymptotic behaviors of the capacities. For particular values of s, p the systems (22) and (28), (31), are dominated by the diagonal terms, for large or small values of κ. This allows the construction of a perturbative scheme and obtains several informations on the asymptotic behaviors of the capacities. The axial symmetry of the problem allows the use of a simplified form of the general Coulombic kernel in (9): ∫ 2π dΦ √ (17) G(r, r0 , x) = 0 r 2 + r02 − 2rr0 cos Φ + x 2 and in terms of this kernel it is immediate to see that the condition of constant potential on the three surfaces can be written in the form ∫ h Lateral surface: dµz0 σ L (z 0 ) [G(a, a, z − z 0 ) ± G(a, a, d + z + z 0 )] −h ∫ a + dµr0 σ B(1) (r0 ) [G(a, r0 , z + h) ± G(a, r0 , z + d − h)] 0 ∫ a + dµr0 σ B(2) (r0 ) [G(a, r0 , h − z) ± G(a, r0 , z + d + h)] = V (18a) 0

∫ Base 1:

h

dµz0 σ L (z 0 ) [G(r, a, z 0 + h) ± G(r, a, d − h + z 0 )] ∫ a + dµr0 σ (1) (r0 ) [G(r, r0 , 0) ± G(r, r0 , d − 2h)] ∫0 a + dµr0 σ (2) (r0 ) [G(r, r0 , 2h) ± G(r, r0 , d)] = V −h

0

∫ Base 2:

h

dµz0 σ L (z 0 ) [G(r, a, h − z 0 ) ± G(r, a, d + h + z 0 )]

−h

(18b)

370

G. Paffuti / Mathematics and Computers in Simulation 166 (2019) 365–381

∫ + 0

∫ + 0

a

dµr0 σ (1) (r0 ) [G(r, r0 , 2h) ± G(r, r0 , d)] a

dµr0 σ (2) (r0 ) [G(r, r0 , 0) ± G(r, r0 , d + 2h)] = V

(18c)

In these equations and in the following we will use the symbols dµr = r dr and dµz = adz for the measures in cylindrical coordinates. The direct projection of these equations on the basis is not convenient due to the bad behavior of the integrand, it is better to transform the kernel (17) with a Fourier transform or with a Hankel transform: ∫ ∞ ∫ ∞ G(r, r0 , x) = 4 dω cos(ωx)I0 (ωr )K 0 (ωr0 ) = 2 dω e−iωx I0 (ωr )K 0 (|ω|r0 ) ; r0 ≥ r (19) 0 −∞ ∫ ∞ G(r, r0 , x) = 2π dω J0 (ωr )J0 (ωr0 )e−ω|x| (20) 0

The first step in the procedure is then to substitute the expansion (14) in (18) and to transform the kernel with (19) where appears σ L and with (20) where appears σ B . The integrals appearing in (18) can be performed using the known properties of Gegenbauer and Jacobi polynomials ∫ a ∫ h J2n−2+ p (ωa) Jn−1+s (ωh) dµr f n (r )J0 (ωr ) = ; . (21) dµz gn (z)eiωz = i n−1 s (ωh) (ωa) p −h 0 These relations hold also with an analytic continuation in ω then can be used also for real exponential functions. At this stage the remaining dependence of Eq. (18a) is on z while the other two depend on r . Projecting the first on gm (z) and the remaining two on f m (r ) with the measures dµz and dµr we have finally the sought linear system ∑

A(0) mn an +

∑

n

∑

n (0) Bmn an +

(1) (1) Bmn bn +

n (0) Cmn an +

n

(2) A(2) mn bn =

∑

(2) (2) Bmn bn =

n

∑

n

∑ n

∑

n

∑

(1) A(1) mn bn +

(1) (1) Cmn bn +

∑

(2) (2) Cmn bn =

n

1 δm1 ≡ ts δm1 2s Γ (1 + s) 1 2pΓ ( p

+ 1)

1 2pΓ ( p

+ 1)

(22a)

δm1 ≡ t p δm1

(22b)

δm1 ≡ t p δm1

(22c)

Here and in the following we will use often for short the notation tx = 1/2x Γ (1 + x). Each matrix in (22) is in correspondence with the analogous term in (18). The explicit calculation is tedious but straightforward, we present the final results in Section 5. We note that in the case of a single conductor this procedure has been developed in [4] and [27]. On the numerical side the generalization given in this paper concerns the case of two conductors, while for a single conductor the only improvement consists in the analytical computation of the matrix elements. Once obtained the solution the charge on the conductor 1 is (∫ h ) ∫ a (1) (2) Q 1 = 2π dµz σ L (z) + dµr (σ B + σ B ) (23) −h

0

Using ∫

h

gn (z)dµz = −h

1 δn1 ≡ ts δn1 ; s 2 Γ (1 + s)

a

∫

f n (r )dµr = 0

1 2pΓ ( p

+ 1)

δn,1 ≡ t p δn1 ;

(24)

one easily finds ( ) 1 Q1 = 2π ts a1 + t p (b1(1) + b1(2) ) (25) a V As explained in the previous section Q 1 /V gives C g1 for equations with both conductors at the same potential and 2C for conductors at opposite potential. Almost all matrices in (22) are given in an explicit analytic form, reported in Section 5. We checked our computations also using a quite accurate version of Boundary Element Method, presented in Section 3.3. In all our computations we take the same number N of polynomials, both for the coordinate z than for the variable r . In usual applications a relatively small N (less than 10) gives accurate results, but as we will be interested in

G. Paffuti / Mathematics and Computers in Simulation 166 (2019) 365–381

371

exploring a quite extreme regime, at very small distances between the electrodes, we used N up to 50 at small values of κ. In the last case the computation can be quite time-consuming and as an help we used a variant of the interpolation method suggested in [19]. At “large” values of κ, surely for κ ≳ 0.1, the results acquire quickly stability with growing N then the error between the computation at small √ N and the stable value can be estimated. We noticed that the errors all lie on a universal curve depending on N k. In the original method, used also in [20] the dependence was on N k. This observation allows, adapting [19], the following extrapolation procedure. (a) We order the values of κ to be computed in decreasing order: κ: κ1 , κ2 . . .. The computation is performed with a maximum number Nmax of polynomials in each variable. Let us call Si (N ) the numerical result obtained for the capacity for the ith term in the above sequence of κ’s using N < Nmax polynomials. (b) For each i the best numerical result is Si (Nmax ). For small κ the extrapolated value Ci for the capacitance is given by [ (√ )] κi Ci = Si (Nmax ) + Ci−1 − Si−1 Nmax (26) κi−1 The necessary values Si can be obtained, if needed, by linear extrapolation. This method is particularly easy to use in the Galerkin method: once the matrix with N = Nmax has been computed, the matrices for lower N are simply extracted by dropping the appropriate lines and columns, and as the computation of the matrix is the time-consuming part of the calculus the advantages are clear. In our case the extrapolation procedure only affects the values of C at very small κ, as will be shown in the next section. In the case of the solution of Love’s equation in [20], an extrapolation similar to (26) was mandatory as N , in that case, was the number of collocation points, and a limit N ≲ 50000 was imposed by memory allocation limits. In Galerkin method the matrices are always relatively small then the memory allocation is not a problem: the extrapolation procedure then is not mandatory but it is useful to save computer-time. In this paper the extrapolated data are mainly used to give an estimation of the error in the computation. Finally we note that for a part of the matrices A(0) , A(2) in (22) we have not been able to find a closed analytical form, then the relative integrals have been computed numerically. The integrand can be highly oscillating and after several experiments we found that the Levin algorithm [10], as implemented in @Mathematica [16], works quite well. 3.2. Single conductor In this case the problem is evidently even (odd) in z then it is more efficient to use a basis with Gegenbauer polynomials of definite parity. For the computation of C1 and D one can use the basis 1 ( z 2 )s−1/2 2s Γ (s)Γ (2n − 1) s ( z ) gn (z) = 1− 2 C (27) 2πah h Γ (2n − 2 + 2s) 2n−2 h The system of equations takes the form ∑ ∑ A(0) A(1) mn an + mn bn = ts δm1 n

∑ n

n (0) Bmn an

+

∑

(1) Bmn bn = t p δm1

(28)

n

The matrices are given in Section 5. The capacity and the quadrupole moment are given by ( ) C1 =2π ts a1 + 2t p b1 a [ ( )2 ( )] C1 D 2π h a1 a2 = −a1 + + + a a 2 2s Γ (1 + s) a 1+s (1 + s)(2 + s) [ ( ( ) ) ] 4π h 2 1 b2 b1 2 − + 2 p Γ (1 + p) a 1+ p (1 + p)(2 + p)

(29a) (29b)

372

G. Paffuti / Mathematics and Computers in Simulation 166 (2019) 365–381

For the computation of the polarizability one proceeds as follows. An external electric field E directed along z produces on the surface of the conductor a potential −E z, the point z = 0 being at the center of the conductor. The surface density must produce the opposite potential, i.e. the right hand side of (9) must be +E z. For a unitary electric field the induced dipole gives directly the polarization, see (11). The charge density is obviously odd in z, then only the odd Gegenbauer polynomials are needed. The basis of odd functions is enumerated by 1 ( z 2 )s−1/2 2s Γ (s)Γ (2n) s ( z ) gn (z) = 1− 2 C (30) 2πah h Γ (2n − 1 + 2s) 2n−1 h while the linear system, projecting the modified equation (9), has the form ∑ ∑ h 1 δm1 A(0) A(1) mn an + mn bn = 1+s 2 Γ (1 + s) a n n (31) ∑ ∑ 1 h (0) (1) Bmn an + Bmn bn = p δm1 2 Γ ( p + 1) a n n The matrices are given in Section 5 and the polarizability is h 1 1 h α = 2π a1 + 4π b1 (32) a3 a 21+s Γ (2 + s) a 2 p Γ ( p + 1) The previous equations (22), (28) and (31) can be specialized to describe particular sub-systems, in particular they can describe two hollow cylinders or one hollow cylinder, described by the matrix A(0) , or the system of two flat discs, where only matrices of the B-type occurs. To assist the reader we list in Section 5 the matrices adapted to these particular cases. 3.3. Boundary element method In the boundary element method (BEM) the surfaces of the conductors are divided into elementary domains (plaquettes) of area Ai and to each plaquette is assigned a charge qi . The distribution of charge must reproduce the assigned potential on the body, i.e. a constant potential for a single conductor. The problem then reduces to a system of linear equations ∑ Vi = Ki j q j (33) j

the sum runs on all plaquettes. Various implementations of the method differ in the choice of the kernel K i j . We choose to consider the charges uniformly distributed on the plaquettes and we perform a mean on the ith plaquette where the potential is computed, i.e. the system (33) takes the form ∫ ∑ 1 ∫ 1 qj (34) Vi = Ai A j x∈Si y∈S j |x − y| j This choice of the kernel is equivalent to a variational calculation with piecewise constant functions, as it is easily seen by considering the energy of the system, see [14]. This implies that the computed capacities are a lower bound of the true result and that the numerical estimates must be growing for finer grids of plaquettes. For axial systems like cylinders we have two kind of surfaces: discs and cylindrical lateral surfaces. The discs are divided in annuli from a starting radius ri to ri + dri , the lateral surface is divided into rings from z i to z i + dz i . Here and in the following the axis z is directed along the symmetry axis of the system. In actual calculations we choose all rings equal, i.e. dz i = dz ∀i, and all annuli with the same area. We have then three different basic forms for the kernel: interaction between two rings of the lateral surface, K (L L) , interaction between two annuli of the bases, K (B B) , and finally the mixed term, K (B L) . In the class of problems under study all single conductors have the same radius a, we can always assume a = 1 for dimensional reasons, all formulas below are given in these units. A straightforward calculation gives [ ] ∫ dz ∫ dz 1 4 4 √ K i(Lj L) = dξ dη K − (35) 2π dz 2 0 (z i + ξ − z j − η)2 (z i − z j + ξ − η)2 0 K is the elliptic integral. The two rings start at coordinates z i , z j .

G. Paffuti / Mathematics and Computers in Simulation 166 (2019) 365–381

For K (B B) we compute: ∫ [ ] 2 1 1 dr j (B B) Ki j = dξ Φ(r j + ξ, ri + dri , z) − Φ(r j + ξ, ri , z) π ai a j 0

373

(36)

where ai = 2ri dri + dri2 . Φ(x, R, z) is proportional to the electrostatic potential of a disc of radius R at a point at distance x from the axis and coordinate z along the symmetry axis: ⎛ ) (√ ⎜ |z| x 2 + z2 + x ⎜ (37) Φ(x, R, z) =2π x ⎜ ⎜ − √ (√ ) ⎝ 2 2 2 2x x +z +x +z ) ( 4Rx ) ( ) ( 4Rx ) R 2 − x 2 − z 2 K (R+x) + (R + x)2 + z 2 E (R+x) 2 +z 2 2 +z 2 √ + π (R + x)2 + z 2 ) ( (√ ) ( (√ )) 2x x 2 +z 2 −x 4Rx | (R+x)2 +z 2 z 2 − (R + x) x 2 + z2 − x Π z2 (

+

√ π (R + x)2 + z 2 )⎞ ( ) ( √ ) ) ( (√ 2x x+ x 2 +z 2 4Rx | (R+x) x 2 + z2 + x + z2 Π − (R + x) ⎟ 2 +z 2 z2 ⎟ ⎟ √ + ⎟ 2 2 π (R + x) + z ⎠ Π, E, K are elliptic integrals, defined by ∫ π/2 ∫ π/2 √ dt √ K (z) = ; E (z) = dt 1 − z sin2 t 0 0 1 − z sin2 t ∫ π/2 dt √ Π(n|m) = 2 0 (1 − n sin t) 1 − m sin2 t The case of planar annuli is given by z = 0. The prefactor 2π x in (37) comes from the integration measure in (36). Finally the factors K (B L) are given by ∫ dz [ ] 1 1 (B L) Ki j = dη Φ(1, ri + dri , z j + η) − Φ(1, ri , z j + η) (38) π dz 2ri dri + dri2 0 Here z j is the distance between the base and the lateral ring, extending from z j to z j + dz. The first argument 1 in Φ is due to our choice of units (a = 1). The integrals in (35), (36) and (38) have to be done numerically: in our computations we used a Gaussian integration routine. Actually the elements K i(Bj B) for z = 0 can be computed analytically but we omit here the long resulting expression as a numerical computation of the integrals gives satisfactory results. Starting from the elements given above it is easy to assemble the whole matrix K i j . For a single conductor one has to solve the system (34) with Vi = 1, ∀i, in this case the sum of charges qi gives directly the capacity. For two conductors it is simpler to break the matrix elements K i j in two parts, Ai j for the self interaction of a conductor, Bi j for the mutual interaction, i.e. B is the part of the matrix K which depends on the distance between the conductors. In this way we realize the decomposition (9) and solving the two systems (Ai j ± Bi j )q j = 1

(39)

374

G. Paffuti / Mathematics and Computers in Simulation 166 (2019) 365–381

we can compute C and C g for the system, as explained in Section 2 and in parallel with the strategy used with the Galerkin approach. In each computation we choose Nr annuli and Nz lateral rings (when both are present). To avoid the introduction of new parameters all computations for full cylinders have been done with a fixed ratio Nr /Nz = 10. As always for BEM the final results have to be extrapolated by intermediate calculations with growing N , we used a quadratic fit in 1/N , here N means Nr for cylinders and discs and Nz for hollow cylinders. The BEM is reasonably accurate and will be a constant check of our calculations with Galerkin method. It has to be noted that for particular systems, like a hollow cylinder, the solution of the problem requires only a number of operations proportional to Nz , in this case the method is very fast. 4. Numerical test In this section we present a numerical test of the method. We consider the case of a couple of parallel flat discs, the only system partially under control from the analytical side. In this case system (22) collapses to a single block ∑ Bi j b j = t p δi1 . (40) j

The matrix Bi j is given by ∫ ∞ ) J2m−2+ p (aω) J2n−2+ p (aω) ( (a) (b) 1 ± e−2hω ≡ Bmn ± Bmn Bmn = 2πa dω p p (aω) (aω) 0

(41)

(a) The matrix B (a) can be computed analytically, with p = 1/2, is diagonal and Bnn = 2π/(4n − 3). The signs ± refer to the computation of C g1 and C respectively and in the two cases:

1 2π t p b1 2 Th analytical expression of (41) in terms of hypergeometric functions is given in (44). The numerical results are presented in Tables 1 and 2. We report for comparison the asymptotic values (5), (6) and the results from a high precision solution [20,22] of Love equation, performed on a grid of 5 × 104 Gaussian points. It is well known that the numerical approximation of this system is quite difficult at low κ, we push our test up to κ = 5 × 10−6 then we somewhat stressed the parameters at our disposal. The boundary elements are computed on a mesh with N B E M = 1200 while the Galerkin method use polynomials of maximum degree N = 500. Up to distances κ ≃ 2×10−4 BEM results are coincident with Galerkin method and Love equation. This provides a quite good cross check for all schemes. At shorter distances κ ≪ 1/N B E M and the method is expected to fail. The first (Galerkin) and third (Love) columns in Tables 1 and 2 are always consistent and comparing the numerical values one appreciates that at very small κ the first column seems more precise, the non extrapolated values are almost identical to the extrapolated values of Love equation. More importantly both calculations reproduce the expected asymptotic behavior (fourth column) within an absolute accuracy of at least 10−8 and a relative precision of 10−12 , for C, i.e. of the order of the expected neglected contributions in (5), (6), i.e. κ 2 log3 κ. We think in effect that the results in Tables 1, 2 are the most accurate results existing in the literature for this system. The computation has clearly also some drawbacks, some of them are related to the extreme regime under study. The calculation of hypergeometric series in this range of parameters proved to be particularly difficult. We used @Mathematica and we need a very high precision, of the order of 500 digits at the smaller values of κ, to avoid rounding errors. We are confident that more expert readers can elaborate more efficient routines. In this particular calculation the computation of the matrix (41) can be fastened by a use of the recursion relations between Bessel functions, but the problem of rounding errors remains. As an intermediate check of the results we computed numerically the integral in (41), in this computation and in general for the integrals of Bessel functions, the Levin algorithm [10] proved to be very useful. C g1 = 2π t p b1 ;

C=

375

G. Paffuti / Mathematics and Computers in Simulation 166 (2019) 365–381 Table 1 Numerical values of C for two discs of radius a = 1. C(B E M) are the results computed with BEM, C(Love) are the results of [20] and the last column the asymptotic expansion at small κ. The lower part of the table gives the extrapolated value obtained with Eq. (26). κ

C/a

C/a (BEM)

C/a (Love)

C/a Eq. (5)

1. 0.7 0.6 0.5 0.4 0.3 0.2 0.15 0.1 0.09 0.08 0.07 0.05 0.04 0.03 0.02 0.01 0.005 0.002 0.001 0.0005 0.0002 0.0001 0.00005 0.00003 0.00002 0.000015 0.00001 0.000005

0.579573860651 0.697478637778 0.762492571650 0.852935151221 0.987494415923 1.20944140950 1.64749629392 2.08070319759 2.93898079847 3.22351516784 3.57841124252 4.03368449989 5.48515774661 6.75077356422 8.85466722892 13.0509956667 25.6031005825 50.6564073534 125.727970864 250.782584053 500.837427208 1250.91012284 2500.96519630 5001.02030718 8334.39426950 12 501.0931901 16 667.7827425 25 001.1483331 50 001.2034705

0.579573860659 0.697478637820 0.762492571714 0.852935151320 0.987494416080 1.20944140976 1.64749629433 2.08070319800 2.93898079787 3.22351516647 3.57841123982 4.03368449476 5.48515772678 6.75077352026 8.85466711422 13.0509952314 25.6030974217 50.6563749532 125.726891645 250.772164433 500.617046770 1263.67616818 48 7316.695495 * * * * * *

0.579573860651 0.697478637778 0.762492571650 0.852935151221 0.987494415923 1.20944140950 1.64749629392 2.08070319759 2.93898079847 3.22351516784 3.57841124252 4.03368449989 5.48515774661 6.75077356422 8.85466722892 13.0509956667 25.6031005825 50.6564073534 125.727970864 250.782584053 500.837427207 1250.91012284 2500.96519631 5001.02030939 8334.39443515 12 501.0947926 16 667.7880979 25 001.1684013 50 001.3060052

0.566663476788 0.689787465066 0.756373847166 0.848283010834 0.984184023190 1.20732225311 1.64638055071 2.08000120515 2.93861918035 3.22321135328 3.57816138551 4.03348452192 5.48504415268 6.75069575380 8.85461962224 13.0509719921 25.6030935132 50.6564052748 125.727970460 250.782583938 500.837427175 1250.91012283 2500.96519630 5001.02030718 8334.39426952 12 501.0931907 16 667.7827444 25 001.1483373 50 001.2034894

k

C-extr.

C (BEM)

C (Love)-extr.

C Eq. (5)

0.0001 0.00005 0.00003 0.00002 0.000015 0.00001 0.000005

2500.96519630 5001.02030718 8334.39426952 12 501.0931907 16 667.7827445 25 001.1483374 50 001.2034893

* * * * * * *

2500.96519631 5001.02030719 8334.39426957 12 501.0931909 16 667.7827448 25 001.1483383 50 001.2034999

2500.96519630 5001.02030718 8334.39426952 12 501.0931907 16 667.7827444 25 001.1483373 50 001.2034894

Finally we make a few comments on the issue of processing time. First we note that in the general case of conductors with non-zero thickness the only comparison term that we can use is the BEM, as the integral Love equation is available only for flat disks. In this case at small κ BEM simply does not work as the fluctuations due to divergent Coulombic forces make the method unstable: in this case we do not know alternatives to the method used in this work. To give an idea of the precision and the times involved we then choose to present in this work the computations on two flat disks, where Love equation gives an efficient alternative. The results of two computations at κ = 0.001, κ = 0.0001 are presented in Table 3. The values for the degree N of the polynomials in Galerkin method, the number of Gaussian points in the numerical solution of Love equation and the number of elements in BEM, Nel , are chosen to give a comparable precision. All procedures are implemented in @Mathematica, running on a 4 cores processor at 4 GHz. We see that both in precision and speed the Galerkin method is better than BEM. The performances are comparable to those of Love equation, which is peculiar for this system. As can be seen in Table 1 Galerkin method is even slightly more accurate at small κ. Let us also note that the computations reported for Love equation in Table 1 are performed with the maximum RAM available on our computer, 64 GB, Galerkin method does not suffer from this type of limitation.

376

G. Paffuti / Mathematics and Computers in Simulation 166 (2019) 365–381 Table 2 Numerical values of C g1 for two discs of radius a = 1. C g1 (B E M) are the results computed with BEM, C g1 (Love) are the results of [20] and the last column the asymptotic expansion at small κ. κ

C g1 /a

C g1 /a (BEM)

C g1 /a (Love)

C g1 /a Eq. (6)

1. 0.7 0.6 0.5 0.4 0.3 0.2 0.15 0.1 0.09 0.08 0.07 0.05 0.04 0.03 0.02 0.01 0.005 0.002 0.001 0.0005 0.0002 0.0001 0.00005 0.00003 0.00002 0.000015 0.00001 0.000005

0.440036084 0.414764706 0.405149049 0.394796318 0.383562858 0.371221845 0.357374396 0.349654371 0.341153441 0.339328286 0.337450709 0.335513814 0.331423372 0.329240989 0.326935911 0.324464479 0.321734448 0.320196629 0.319157117 0.318768559 0.318556765 0.318417917 0.318367413 0.318340405 0.318328974 0.318323022 0.318319957 0.318316805 0.318313520

0.440036084 0.414764706 0.405149049 0.394796318 0.383562858 0.371221844 0.357374396 0.349654371 0.341153441 0.339328286 0.337450710 0.335513814 0.331423372 0.329240990 0.326935912 0.324464481 0.321734457 0.320196684 0.319157907 0.318771090 0.318559937 0.318413929 0.318340041 * * * * * *

0.440036084 0.414764706 0.405149049 0.394796318 0.383562858 0.371221845 0.357374396 0.349654371 0.341153441 0.339328286 0.337450710 0.335513814 0.331423372 0.329240989 0.326935911 0.324464479 0.321734448 0.320196629 0.319157117 0.318768559 0.318556765 0.318417917 0.318367413 0.318340405 0.318328974 0.318323022 0.318319957 0.318316805 0.318313521

0.426963171 0.407015741 0.399029094 0.390194152 0.380339133 0.369204064 0.356347559 0.349024262 0.340840247 0.339067597 0.337238542 0.335345994 0.331330829 0.329178824 0.326898814 0.324446660 0.321729426 0.320195232 0.319156864 0.318768490 0.318556746 0.318417914 0.318367412 0.318340405 0.318328974 0.318323022 0.318319957 0.318316805 0.318313521

Table 3 Results at intermediate precision for C in a two disk capacitor and their processing times. The “exact” results are those of Table 1. BEM is clearly unsuitable for the lower κ. κ = 0.001 Galerkin Love eq. BEM

N = 100 N Gauss = 104 Nel = 500

κ = 0.0001

C/a

Time

C/a

Time

250.782584052 250.782584053 250.739456589

16 s 10.2 s 64.7 s

2500.96512789 2500.96679700 −17 221.38

14.5 s 10.2 s 64.7 s

Table 4 Results at intermediate precision for C and C g1 in a two disk capacitor and their processing times for κ = 0.5, τ = 0.05. Galerkin Galerkin BEM BEM

N = 10 N =5 Nel = 500 Nel = 200

C/a

C g1 /a

Time

0.8982849 0.8982832 0.897529 0.897136

0.4148418 0.4148412 0.414767 0.414717

24.7 s 6.4 s 270. s 44. s

In the case of two thick conductors the comparison is almost impossible as BEM loses quickly in precision for small κ. In Table 4 we give an example for κ = 0.5 and τ = 0.05. In comparing the results we have to add that for BEM it is necessary to perform the computation at several N and perform an extrapolation, this requires an additional computer time.

377

G. Paffuti / Mathematics and Computers in Simulation 166 (2019) 365–381

5. Details on the matrix elements 5.1. Two conductors The cylinders have radius a and length 2h. d is the distance between the centers, then the distance between the nearest surfaces is d − 2h = ℓ ≡ a κ. In actual computations clearly one can put a = 1. The signs ± apply to the computation of C g or C respectively, as explained in previous sections. A(0) mn with

∫ =a

∞

Jm−1+s (ωh) Jn−1+s (ωh) m+n−2 dω 4I0 (ω)K 0 (ω) i · X mn (ωh)s (ωh)s 0 { ±i sin(ωd) m + n odd X mn = (−1)m−1 ± cos(ωd) m + n even

m−1 A(1) a mn = 2π (−1)

∫

(2) Bmn (2) Cmn

(42b)

] Im−1+s (ωh) J2n−2+ p (ωa) [ −ωh e ± (−1)m−1 e−ω(d+h) p (ωa) (ωh)s

(42c)

A(2) nm

dω J0 (ωa) ∫

] Im−1+s (ωh) J2n−2+ p (ωa) [ −ωh e ± e−ω(d−h) p (ωa) (ωh)s

dω J0 (ωa)

∞

0

(1) Bmn

∞

0 (0) Cmn =

(0) Bmn = A(1) nm ;

A(2) mn = 2πa

∫

∞

] J2m−2+ p (ωa) J2n−2+ p (ωa) [ 1 ± e−ω|d−2h| p p (ωa) (ωa) ∫0 ∞ ] J2m−2+ p (ωa) J2n−2+ p (ωa) [ −2ωh dω = 2πa e ± e−ω|d| ; p p (ωa) (ωa) ∫0 ∞ ] J2m−2+ p (ωa) J2n−2+ p (ωa) [ dω = 2πa 1 ± e−ω(d+2h) p p (ωa) (ωa) 0

= 2πa

(42a)

(42d)

dω

(1) (2) Cmn = Bmn

(42e) (42f)

We see that each matrix is composed by two parts, the second one depends on d while the first one is the only which will be needed for the computation of the capacities of a single conductor (see below). Moreover in a computation performed on several distances this part of the matrix can be computed only once. Accordingly we will denote the separate parts of the first matrix with the notation A(0,a) , A(0,b) and similarly for the other cases. All the matrices elements of type (a) can be computed analytically for generic values of the parameters s, p of Gegenbauer end Jacobi polynomials. Their expression is the following, we put a = 1 to simplify the formulas: 1 · A(0,a) mn = ( π1 h 1 ) 1 2 , 2 (−m − n + 3), − |m−n| + s + 21 , 12 (m + n + 4s − 1), 21 (|m − n| + 2s + 1) 3,3 2 G 5,5 | 0, 0, s, 0, s + 12 h2 2π m−1 (1,a) A(1,a) (−1)m−1 F J J I E(m − 1 + s, 2n − 2 + p, p + s, h) ; A(2,a) Amn mn = mn = (−1) hs (1,a) (2,a) = 2π F J J E(2m − 2 + p, 2n − 2 + p, 2 p, 2h) Bmn = 2π F J J (2m − 2 + p, 2n − 2 + p, 2 p) ; Bmn (2,a) Cmn = 2π F J J (2m − 2 + p, 2n − 2 + p, 2 p) (1,b) Bmn = 2π F J J E(2m − 2 + p, 2n − 2 + p, 2 p, d − 2h) (2,b) Bmn = 2π F J J E(2m − 2 + p, 2n − 2 + p, 2 p, d) (2,b) Cmn = 2π F J J E(2m − 2 + p, 2n − 2 + p, 2 p, d + 2h)

Here G is the Meijer G-function, the functions F J J, F J J E, F J J I E denote respectively integrals of two Bessel functions and a power, two Bessel functions an exponential and a power, and finally two Bessel functions a Bessel I function, an exponential and a power. Their explicit general expression is given below, after the presentation of

378

G. Paffuti / Mathematics and Computers in Simulation 166 (2019) 365–381

the matrices. We have used the notations of reference [5] and we used the Mathematica software [16] to perform some of the algebraic computations. Previous matrices can be specialized to particular systems. Two hollow cylinders - In this case only the matrix A(0) mn above, (42a), is needed. The natural parameter for Gegenbauer polynomials is s = 0 in this case. Two flat discs - Only a combination of matrices B, C above appear. We give here the explicit numerical and analytical expression (with a = 1 in the latter): ∫ ∞ ) J2m−2+ p (aω) J2n−2+ p (aω) ( (a) (b) 1 ± e−2hω ≡ Bmn ± Bmn (43a) Bmn = 2πa dω p p (aω) (aω) 0 { } Bmn = 2π F J J (2m − 2 + p, 2n − 2 + p, 2 p) ± F J J E(2m − 2 + p, 2n − 2 + p, 2 p, 2h) (43b) We remember that, in this case, with a = 1, d = 2h is the distance between the discs. The natural choice of the parameter p is p = 1/2. In this case the first part of the matrix, the one non depending on h, is diagonal and (a) Bnn = 2π /(4n − 3). The system to be solved and the capacities are given by ∑

Bi j b j = t p δi1 ;

C g1 = 2π t p b1 ;

j

C=

1 2πt p b1 2

(44)

depending on whether one chooses the plus or minus sign in (43).

5.2. A single conductor In the case of a single cylinder of radius a and length 2h it is convenient to write separately the contribution of even and odd Gegenbauer polynomials, in the first case we can compute the capacity and the quadrupole moment (for unit charge). In the second case we can compute the polarizability. According to our previous notations the length of the cylinder is L = 2h. In the following all matrix elements are given for cylinder of unit radius. For the computation of C and D the system to be solved is given in (28), where: [ ] ∫ ∞ J2m−2+s (ωh) J2n−2+s (ωh) n+m dω (−1) A(0) = a 4I (ωa)K (ωa) (45a) 0 0 mn (ωh)s (ωh)s 0 A(1) mn = a

∫

∞

dω 4π e−ωh J0 (ωa)

0

(1) Bmn

∫ =a 0

∞

dω 2π (1 + e−2ωh )

I2m−2+s (ωh) J2n−2+ p (ωh) (0) = 2Bnm (ωh)s (ωh) p

(45b)

J2m−2+ p (ωa) J2n−2+ p (ωa) (ωa) p (ωa) p

(45c)

The prefactor 2 in B (0) is due to our normalization to the charge of a single basis of the cylinder. The analytical expression can obtained from matrix of type (a) essentially by the replacement 2m → 2m − 1 in the indices relative to Gegenbauer polynomials, in any case the explicit expression is ⎛ ⎞ 1 5 1 3 1 , −m − n + , − |m − n| + s + , m + n + 2s − , |m − n| + s + ⎟ 1 3,3 ⎜ 1 2 2 2 2 2⎠ A(0) · G 5,5 ⎝ 2| mn = 1 πh h 0, 0, s, 0, s + 2 4π (1) Amn = s F J J I E(2m − 2 + s, 2n − 2 + p, p + s, h) h { } (1) Bmn = 2π F J J (2m − 2 + p, 2n − 2 + p, 2 p) + F J J E(2m − 2 + p, 2n − 2 + p, 2 p, 2h)

379

G. Paffuti / Mathematics and Computers in Simulation 166 (2019) 365–381

For the computation of α the system to be solved is given in (31), where: [ ] ∫ ∞ J2m−1+s (ωh) J2n−1+s (ωh) (0) n+m Amn = a dω (−1) 4I0 (ωa)K 0 (ωa) (ωh)s (ωh)s ∫0 ∞ J (ωh) I (ωh) 2n−2+ p 2m−1+s (0) A(1) dω 4π e−ωh J0 (ωa) = 2Bnm mn = a (ωh)s (ωh) p 0 ∫ ∞ J2m−2+ p (ωa) J2n−2+ p (ωa) (1) Bmn = a dω 2π (1 − e−2ωh ) (ωa) p (ωa) p 0

(46a) (46b) (46c)

The analytical results for the integrals are: ⎛ 3 1 1 1 1 3,3 ⎜ 1 2 , −m − n + 2 , − |m − n| + s + 2 , m + n + 2s − 2 , |m − n| + s + (0) Amn = G ⎝ | 1 π h 5,5 h 2 0, 0, s, 0, s + 2

⎞ 1 2⎟ ⎠

4π F J J I E(2m − 1 + s, 2n − 2 + p, p + s, h) hs { } = 2π F J J (2m − 2 + p, 2n − 2 + p, 2 p) − F J J E(2m − 2 + p, 2n − 2 + p, 2 p, 2h)

A(1) = B (1)

As in the previous section the matrices can be specialized to describe a hollow cylinder. In both cases only the matrices A(0) mn has to be used. 5.3. Integrals

F J J I E(µ, ν, λ, α) =

∫

∞

dt J0 (t)Iµ (αt)Jν (t)t −λ e−αt

(47)

0

{ ( ) ( ) 16αΓ λ + 12 Γ 41 (−2λ + 2ν + 1) 2−λ−5 = √ 3/2 ) ( )2 · ( πα Γ 14 (2λ − 2ν + 3) Γ 41 (2λ + 2ν + 3) ( λ 1 λ 3 1 µ 3 µ µ 1 µ 3 + , + , − , − , + , + ; 6 F5 2 4 2 4 4 2 4 2 2 4 2 4 1 λ ν 3 λ ν 3 λ ν 3 λ ν 3 1 , − + , − + , + + , + + ;− 2 2 2 2 4 2 2 4 2 2 4 2 2 4 α ( ) ( ) (2µ − 1)(2µ + 1)Γ λ + 23 Γ 14 (−2λ + 2ν − 1) − · ( ) ( )2 Γ 14 (2λ − 2ν + 5) Γ 41 (2λ + 2ν + 5) ( λ 3 λ 5 3 µ 5 µ µ 3 µ 5 F + , + , − , − , + , + ; 6 5 2 4 2 4 4 2 4 2 2 4 2 4

)

) 5 λ ν 5 λ ν 5 λ ν 5 1 3 λ ν , − + , − + , + + , + + ;− 2 2 2 2 4 2 2 4 2 2 4 2 2 4 α ( ) 1 4λ−ν+2 α λ−ν+ 2 Γ λ − ν − 12 Γ (−λ + µ + ν + 1) + · Γ (ν + 1)Γ (λ + µ − ν) ( ν 1 ν λ µ ν 1 λ µ ν λ µ ν 1 λ µ ν F + , + 1, − − + + , − − + + 1, − + + + , − + + + 1; 6 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 )} λ ν 3 λ ν 5 1 1, − + + , − + + , ν + 1, ν + 1; − 2 2 2 4 2 2 4 α

380

G. Paffuti / Mathematics and Computers in Simulation 166 (2019) 365–381

F J J (µ, ν, λ) =

∫

∞

dt Jµ (t)Jν (t)t −λ

0

( ) 2−λ Γ (λ)Γ 12 (−λ + µ + ν + 1) ) ( ) ( ) = (1 Γ 2 (λ + µ − ν + 1) Γ 21 (λ − µ + ν + 1) Γ 12 (λ + µ + ν + 1) F J J E(µ, ν, λ, α) =

(48)

∞

∫

dt Jµ (t)Jν (t)t −λ e−αt = 2−µ−ν Γ (µ + ν + 1)α λ−µ−ν−1 Γ (−λ + µ + ν + 1)

0

˜ 4 F3

(

1 1 1 1 (µ + ν + 1), (µ + ν + 2), (−λ + µ + ν + 1), (−λ + µ + ν + 2); 2 2 2 2 ) 4 µ + 1, ν + 1, µ + ν + 1; − 2 α

(49)

Here p Fq is the generalized hypergeometric function, and F˜ its regularized form. 6. Conclusion In this work we present a method based on Galerkin expansion for the computation of the capacity matrices for circular thick electrodes. This method can be used for the estimation of the forces between the conductors, a problem of obvious importance in the field of Micro Electro-Mechanical Systems (MEMS). A parallel construction of a boundary elements procedure, with an integration of the Coulomb kernel both on the source and field patches provide an alternative tool and a check of the more accurate Galerkin procedure. We test the method on the two flat discs capacitor and obtain a very good agreement with the known asymptotic formulas for this system. References [1] C. Atkinson, F. Leppington, The asymptotic solution of some integral equations, IMA J. Appl. Math. 31 (3) (1983) 169–182. [2] G. Carlson, B. Illman, The circular disk parallel plate capacitor, Amer. J. Phys. 62 (12) (1994) 1099–1105. [3] W. Chew, J. Kong, Microstrip capacitance for a circular disk through matched asymptotic expansions, SIAM J. Appl. Math. 42 (1982) 302–317. [4] S. Falco, G. Panariello, F. Schettino, L. Verolino, Capacitance of a finite cylinder, Electr. Eng. 85 (4) (2003) 177–182. [5] I.S. Gradshteyn, I.M. Ryzhik, Table of integrals, series, and products. edited by a. jeffrey and d. zwillinger, 2007. [6] V. Hutson, The circular plate condenser at small separations, in: Mathematical Proceedings of the Cambridge Philosophical Society, Vol. 59, Cambridge University Press, 1963, pp. 211–224. [7] T. Itoh, R. Mittra, A new method for calculating the capacitance of a circular disk for microwave integrated circuits (short papers), IEEE Trans. Microw. Theory Tech. 21 (6) (1973) 431–432. [8] G. Kirchhoff, Zur theorie des condensators, in: Monatsber. Der Akad. Der Wiss. Zu Berlin, 1877, pp. 144–162. [9] L. Landau, E. Lifshitz, Electrodynamics of Continuous Media, Pergamon Press, Oxford, 1984. [10] D. Levin, Analysis of a collocation method for integrating rapidly oscillatory functions, J. Comput. Appl. Math. 78 (1) (1997) 131–138. [11] R. Love, The electrostatic field of two equal circular co-axial conducting disks, Quart. J. Mech. Appl. Math. 2 (1949) 428–451. [12] F. Maccarrone, G. Paffuti, Capacitance and potential coefficients at large distances, J. Electrost. 83 (2016) 22–27. [13] F. Maccarrone, G. Paffuti, Some comments on the electrostatic forces between circular electrodes, J. Electrost. 84 (2016) 135–142. [14] F. Maccarrone, G. Paffuti, Capacitance and forces for two square electrodes, J. Electrost. 89 (2017) 20–29. [15] F. Maccarrone, G. Paffuti, Capacitance and forces for thick circular electrodes, J. Electrost. 94 (2018) 30–37. [16] Mathematica 11.0, Wolfram Research Inc., Champaign, Illinois, 2016. [17] J. Meixner, The behavior of electromagnetic fields at edges, IEEE Trans. Antennas and Propagation 20 (4) (1972) 442–446. [18] J. Nicholson, Oblate spheroidal harmonics and their applications, Phil. Trans. Roy. Soc. Lond. Ser. A 224 (1924) 49–93. [19] M. Norgren, B. Jonsson, The capacitance of the circular parallel capacitor obtained by solving the love integral equation using an analytic expansion of the kernel, Prog. Electromagn. Res., PIER 97 (2009) 357–372, 209. [20] G. Paffuti, Numerical and analytical results for the two discs capacitor problem, in: Proc. R. Soc. A, Vol. 473, The Royal Society, 2017, p. 20160792. [21] G. Paffuti, Results for Capacitances and Forces in cylindrical systems, ArXiv e-prints, Jan. 2018. [22] G. Paffuti, E. Cataldo, A. Di Lieto, F. Maccarrone, Circular plate capacitor with different discs, in: Proc. R. Soc. A, Vol. 472, The Royal Society, 2016, p. 20160574. [23] A.H. Scott, H. Curtis, Edge correction in the determination of dielectric constant, J. Res. 22 (1939) 747–775.

G. Paffuti / Mathematics and Computers in Simulation 166 (2019) 365–381 [24] [25] [26] [27] [28]

S.J. Shaw, Circular-disk viscometer and related electrostatic problems, Phys. Fluids 13 (1970) 1935–1947. P.P. Silvester, R.L. Ferrari, Finite Elements for Electrical Engineers, Cambridge university press, 1996. I.N. Sneddon, Mixed Boundary Value Problems in Potential Theory, North-Holland, 1966. L. Verolino, Capacitance of a hollow cylinder, Electr. Eng. 78 (4) (1995) 201–207. L. Wigglesworth, Comments on ‘circular-disk viscometer and related electrostatic problems’, Phys. Fluids 15 (1972) 718.

381