Perturbative Approach for the Analysis of Charge Distribution on Arbitrarily Shaped Conductors

In this paper, we analyze the electrostatic charge distribution on arbitrarily shaped conductor surfaces. Following a perturbative approach, we derive an approximate analytical formulation of the problem. We start from the known case of a conducting ellipsoid, we adopt a deformed ellipsoidal coordinate system, and we search for the zero-order approximated solution of the problem. We also focus on arbitrary-shaped thin foils, showing that the charge density is divergent on their borders. We then define the applicability range of the proposed approach expressing the contour equation as the Fourier series. Finally, we present a detailed error analysis for several polygonal contours, comparing the analytical results with those obtained via a numerical analysis based on the Finite Element Methods (FEM).


I. INTRODUCTION
The problem of finding the electrostatic charge distribution on a charged conductor is a fascinating problem that many classical electrodynamics textbooks solve in several analytical cases [1]- [3]. Still today it is attracting the attention of the researchers. In [4], [5] the authors study the relationships between the electrostatic field near the conducting surface and curvature radius of the surface itself. In [6], [7] the authors perform a numerical analysis of electrostatic problems. In [8], [9] the authors, respectively, study analytically and numerically the capacitance of different arbitrarily shaped conductors. Several practical problems are related to the identification of the charge distribution and possibly its optimization, on thin surface domains, for instance, in energy harvesting and designing devices [10], [11], semiconductors [12], photovoltaics, and optoelectronics [13], [14]. In electrical accumulators, such as lithium batteries, the study of the charge density distribution has a crucial role in developing new techniques to reduce dendrite growth [15], [16]. Many devices used in electrical and electronic engineering are sketched with an ideal geometrical The associate editor coordinating the review of this manuscript and approving it for publication was Jenny Mahoney .
shape. The knowledge of the charge distribution of a generic geometry of the conductor can allow for a better designing the conductor geometries and for reducing the possibility of occurrence of electrostatic discharge [17]- [19] A typical example is a capacitor considered as two infinite metallic foils separated by a distance d. This idealization of a capacitor does not take into account the intense electric field generated near the edges of the conducting surfaces. The problem of electrostatic forces between two finite-size conductors of different geometries has been investigated in several works [20]- [24]. Recently, in [25], the authors consider a battery with finite-size electrodes, making the conjecture that the border effect can be relevant in the dendrite formation. Being of finite size and assumed as two twodimensional disks, the electrodes have a not uniform charge density taking infinite values at the border of each disk. More generally, any device has a finite size, and the knowledge of the charge density distribution is crucial. The proposed model has been supported by experiments [26], so confirming that the density charge distribution on the battery electrodes plays an important role in dendrite formation. While uniformly charged two-dimensional shapes are studied in the literature, iso-potential surfaces with non-uniform charge distribution are a hard task, and closed formulas for the distribution are VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ known in a few cases [1]- [3]. In [27]- [30] the authors present interesting methods based on series development. The series method may be difficult to apply when in the solutions, there are divergent quantities. Typically, to reproduce divergent quantities, an increasing number of series terms is necessary. Similar difficulties can arise in numerical approaches such as, for example, the Finite Element Method (FEM). It is an intrinsic limitation for a numerical simulation to reproduce a divergence. The approach that we propose can overcome the difficulty of dealing with divergent quantities since the leading term, zero-order term, in thin foils, is divergent on the foil border. In the present work, we propose a general approach to the problem of an analytical expression of the potential, the electrical field, and the electrical charge density of a metallic surface with an arbitrary shape. In particular, we will focus on thin foil surfaces. The paper is organized as follows: In Sec. II, we will show the analytical approach to the full problem, i.e., the evaluation of the electric potential and the charge density on an arbitrarily shaped conductor. In Sec. III, we will apply the results of the previous section to thin foils. In Sec. IV, we will show the applicability of our approach. In Sec. V, we will perform a numerical analysis and check on the analytical formulas for thin foils with polygonal shapes, comparing the results with the simulations performed by the software Ansys. In Sec. VI we will discuss the numerical results. Finally, in Sec. VII, we will summarize the paper results.

II. ELECTRICAL POTENTIAL GENERATED BY A CONDUCTOR SURFACE
Our goal is to find an analytical (approximate) solution of where S is the conductor surface and φ 0 is the value of the potential at the surface. Since we are considering an arbitrary shape for the conductor surface, we need to adopt a convenient set of coordinates, For this purpose, we consider the coordinate system generated by the following equation [2] r 2 a 2 (1 + εf (θ )) 2 + u + where r 2 = x 2 + y 2 . It is straightforward to check that u = 0 describes the deformed ellipsoid where its section has an arbitrary shape given by a(1 + εf (θ )) (see the example of Fig. 1). Keeping in mind thin foils, we may focus on the confocal oblate spheroids coordinates, i.e., a(1 + εf (θ)) > c so we may take the limit c → 0 and end up into the twodimensional thin foil case. Solving for u we have where R 2 = r 2 + z 2 , a(θ ) = a(1 + εf (θ )) and We can choose as third coordinate the polar angle θ . Without ambiguity from now on we set u = u 1 and u 2 = v. The inverse relationships are To evaluate the laplacian of the potential in the u, v, θ coordinates, we may use the general formula given by [31] where is the determinant of the metric tensor, g ik , and g ik are the inverse matrix elements of g ik . Also, according to Einstein's summation convention, repeated suffix implies summation. To find the metric tensor associated with these coordinates, we write the infinitesimal distance between two points with g ik = g ki , u 1 = u, u 2 = v, u 3 = θ . Plugging Eqs. (4) and (5) into Eqs. (7) and (8). We obtain the metric tensor elements writing the square of the infinitesimal length of a curve in the u, v, θ coordinates where for the sake of compactness, we did not explicitly write the off-diagonal terms, δg ik . From Eq. (10), we deduce the elements of the metric tensor From here on, the parameter ε will be considered a perturbative parameter, i.e., ε 1. To zero ε order, we seek a solution that is a function of the coordinate u, namely the coordinate that vanishes at the surface of the conductor. The Laplace equation for the electric potential φ(u, v, θ) is given by While the coordinates u, v do not have crossed terms, g uv = 0, in principle we should consider the fact that g θ k = 0. But the contribution of these terms and the θ derivative is of the ε order (see Eqs. (14) and (15)). We also stress that we cannot neglect ε in the function a(θ ) in the denominator of h 1 because the variable u can take the zero value since u = 0 is the surface of the conductor. Therefore, in the denominator, the term containing ε cannot be neglected. Solving Eq. (16) with the boundary condition φ → 0 for u → ∞ and φ = constant for u = 0, we have and the charge density on the surface of the conductor, to zero-order, is where 0 is the vacuum permittivity.

III. CHARGE DENSITY ON THIN FOILS
In this section, we will focus now on the case of a thin foil with an arbitrary border. The foil's electric potential is obtained taking the limit c → 0 of Eq. (17). The surface density is obtained taking the limits c → 0, u → 0 and z → 0 of Eq. (18). Adopting the more visualizable cylindrical coordinate system r, θ, z, for the potential we have where For the density we have We infer that at the border of the conductor, r = a(θ ), the density is divergent. When ε = 0 then a(θ ) = a[+εf (θ )] = a and we recover the exact result of a disk of radius a. As we will see in Sec. V, for ε 1 the zero approximation is already satisfactory.
We now evaluate the first correction to the potential, and consequently, to the charge density. To further simplify our calculation, we use cylindrical coordinates. We focus on solving the Laplace equation near the surface, i.e., at z = 0, so finding the correction to the charge density. We may write where φ 1 is of ε order. The potential φ 0 satisfy the Laplace equation in the r, z variables while the derivative respect to θ generates a ε order term. Developing φ 0 and φ 1 near z = 0 we have We stress that the problem is symmetric with respect to the z coordinate, and the potential is an odd function of z. In principle, we should also consider the term proportional to z 3 but, for the moment, we assume that it is negligible with respect φ 1 . Assuming a Taylor expansion in r for the function (r, θ), we write and plugging it into Eq. (23) we have ∞ n=0 4n 2 r 2n y n (θ) + r 2n ∂ 2 ∂θ 2 y n (θ ) where we defined h n (θ ) ≡ a(θ ) −2n−1 . Solving for n = 0 we have while for n = 0 we have where h (θ ) means second derivative with respect to θ and the integration constants have been selected in such way that y n (θ ) → 0 for ε → 0. After long but straightforward algebra we obtain For sake of simplicity, we assume f (0) = 0 which, in turn, implies h (0) = 0. The condition f (0) = 0 is satisfied by all polygons with an even number of sides. Performing the sum we have where U (x) is the step function. Remembering that the correction to the potential near z = 0 is φ 1 = z n r 2n y n (θ) we obtain for the density correction By direct inspection of Eq. (30), we may verify that for ε → 0 then σ (1) (r, θ) → 0. To take into account the term proportional to z 3 , previously neglected, we introduce the effective constant B, instead of the constant A multiplying σ (1) (r, θ). Once the density correction has been determined, we can write for the electric potential φ(x, y, z) where R 2 = x 2 + y 2 + z 2 , the extra factor 2 takes into account the contribution of the upper and lower conductor surface, and with where, as previously stated, the effective constant B multiplies σ (1) (r, θ) instead of the constant A. Also A and B have been redefined in such a way to include the constant 2 4π 0 . By construction, φ(x, y, z) given by Eq. (31) satisfies the Laplace equation and ∂ z φ(x, y, z) | z=0 gives the charge density at the surface of the conductor, i. e. σ (r, θ) = σ (0) (r, θ)+ σ (1) (r, θ). In Sec. V we will numerically check that φ(x, y, z), at z = 0 and in the region r ≤ a(θ), is, with good approximation, constant.

IV. APPLICABILITY OF THE APPROXIMATION
So far, we have not given details on the applicability of the approximation in the presence of thin foils with arbitrary border. It is important to note that, despite a shape that visually appears relatively far from a circular shape, what is relevant in our calculation is the first non-constant term of a(θ ) = a(1 + εf (θ )), i.e. εf (θ ). If ε 1, then our approach applies. In general, the border is described by a function where not necessarily identifying the ε parameter is straightforward. We can give a criterion to identify the parameter using the Fourier series of the border. To do that, we write where c * k is the complex conjugate of c k . Let us definek the index of the first Fourier coefficient that does not vanish.
Then, as a first approximation, we may write a(θ) as the sum of the first two non-vanishing terms where we use the polar representation of the Fourier coefficients ck = |c|k exp[iφk ]. Our approach applies if 2|c|k /c 0 1. Indeed, being c k the coefficient of a convergent series, then in general we have that c k < ck for k >k.
To better clarify this point, let us consider the case of a square. Its border is described by that represents a square of side a = √ 2. The first nonvanishing terms are given by c 0 ≈ 0.79, c 4 = c * 4 ≈ −0.055 and we can write Consequently we may set ε = 2 c 4 /c 0 ≈ 0.14. We deduce that a square can be considered a deformed circle, even if the first-order perturbative approach is not accurate. We will return to the square conductor study in the next section to perform numerical analysis on surfaces with a polygonal border.

V. NUMERICAL ANALYSIS
In this section, we will check our result by performing a numerical analysis on a thin foil with a polygonal edge with n sides. The function that describes a n-sided polygon is The main goal of this section is the numerical evaluation of the percent error defined as E = 100(φ − φ S )/φ S at z = 0 in the region r ≤ a(θ ) where the potential φ takes the constant value φ S . Without loss of generality we may set φ S = 1. In order to perform the numerical analysis, the coefficients A and B of Eqs. (32) and (33) need to be evaluated. We calculate the function φ in a D-elements computational domain (x i , y i ) belonging to the polygon. We define M 1 (x i , y i ) and M 2 (x i , y i ) as where r 2 i ≡ x 2 i + y 2 i and A and B are determined by solving the system . . . . . . (41) The linear system of Eq. (41) is over-determined, and A and B are calculated by minimizing the 2-norm N (2) given by To evaluate M 1 (x i , y i ) and M 2 (x i , y i ) we perform a numerical integration by means of Matlab R [32]- [34]. First of all, we generate a two-dimensional numerical domain, according to Eq. (38). Since the numerical values of the integrals defining M 1 and M 2 are divergent along the border of the domains, the distribution of the couples (x i , y i ) inside the polygonal domain is crucial for the evaluation of A and B. Therefore, we randomly selected 6×10 3 couples for each case to provide a uniform point distribution and a minor bias in the A and B estimation. Also, we adopted a minimum distance between the points to increase the uniformity of their distribution. The adopted number of coordinates, at which M 1 and M 2 are evaluated, is a balance between the accuracy calculating A and B and the computational resources availability. In Fig. 2 we show, as example, the domain for n = 4 (i.e. square domain) resulting from the aforementioned discretization. Being the surface is an electrical conductor, the domain r ≤ a(θ ) is expected to be at the same potential φ 0 = 1, and, ideally, A and B are expected to satisfy AM 1 (x i , y i ) + BM 2 (x i , y i ) − 1 = 0 ∀i. In the practical case, the analytical approximation and the numerical errors, mainly due to the adopted integral algorithm tolerance [33], [34], lead to different values from the expected one. With the choice φ S = 1, the percent error E for each domain position is given by The A and B coefficients are calculated for a polygon with n sides and n ranging in n ∈ [4,14] and minimizing N (2) given by Eq. (42). Average, median, maximum positive VOLUME 9, 2021 and negative errors based on E are further compared for all the tested cases. Finally, we compare the analytical results with the numerical simulation performed using the software Ansys. More precisely, we compare the numerical simulation with the zero-order of the perturbative approach given by where the surface potential φ S , for numerical reasons, has been set 100 volts instead of the unit. The thickness of the thin octagonal foil is 0.01 mm, and the radius of the circumscribed circumference, i.e., the distance center-vertex of the octagon, is 1 mm. The agreement between the analytical formula (44) and the numerical analysis is very good, as shown in Figs. 3 and 4. It is worthy of stressing that, being the charge density divergent on the octagon border, the agreement between analytical formulation and numerical calculations ceases to hold in the vicinity of the border. Red line the numerical evaluation. We compared the results along the line θ = 0, i.e., the x axis. For numerical reason the potential on the surface φ S has been set 100 volts. The thickness of the foil is 0.01 mm, and the radius of the circumscribed circumference is 1 mm.

FIGURE 4.
Charge density: Blu line the analytical result given by (32 ). Red line the numerical evaluation. We compared the results along the line θ = π/8, i.e., from the origin to the vertex of the octagon. For numerical reason the potential on the surface φ S has been set 100 volts. The thickness of the foil is 0.01 mm, and the radius of the circumscribed circumference is 1 mm.

VI. DISCUSSION OF THE NUMERICAL RESULTS
The potential φ resulting from the numerical integration of Eq. (31), with the A and B coefficient obtained according to the aforementioned method, is calculated, and the error E derived from Eq. (43) is shown in Fig. 5 for different geometries (φ S = 1). In particular, it is shown that the value of φ is underestimated with respect to the attended value φ S = 1 in the corner region in all the cases. At the same time, the error E practically vanishes in the center of the domain. The error E is significantly reduced with the increase of the polygon sides, as shown in Figs. 5 a-d. This behavior is expected since, increasing n, the influence of the corners is reduced. The domain tends to a circular shape and Eq. (19) with a(θ ) = a gives the exact expression for potential φ(r, θ). A quantitative analysis is performed on the maximum positive and negative percent error E, the median and mean error along with the numerical domain, which is shown in Fig. 6.  The error E has a descent trend as a function of n, while the domain tends to the circular one. By comparing Figs. 5 and 6, we observe that the largest error is located on the corner area, where it reaches value −18% on the square domain case. A slight overestimation of φ occurs on the side area within a 9% range. Regarding the numerical calculations, the largest error amount is found on the domain borders, where Eqs. (32) and (33 ) diverge, thus introducing an underestimation of φ during the numerical integration. The presented calculations show an average underestimation error within 0.2% range and a median error within 0.35% range for all the tested cases, thus confirming the capability of the current method to perform a satisfactory evaluation of the electrical potential φ and the charge density σ on the chosen domains. The obtained values of A, B, and the norm N (2) for the different domains are shown in Tab. 1.  (2) for the polygonal cases. The value of A is of the order of π −2 ≈ 0.1013.  According to Eq. (31), at x = y = z = 0, the integral associated to the σ (0) term takes the value Aπ 2 for any border a(θ ). The value of M 1 is expected approximatively π 2 , and consequently the value of A is of the order of π −2 ≈ 0.1013 as consistently shown in Tab. 1. The coefficient B decreases as the number of sides increases, which is also expected. Indeed, while n increases, the domain tends to a circle, reducing the corrective term's contribution to zero. However, it is worthy to notice that the low contribution of the M 2 (x i , y i ) term to the potential φ leads to higher sensitivity in B calculation through the solution of the over-determined system shown in Eq. (41). Finally, Figs. 7b and 7d show the charge density behavior on a sample line i.e. y = 0 in order to highlight its divergent behavior on the borders.

VII. CONCLUSION
In this paper, we presented an analytical formulation of the problem of electrostatic charge distributions on arbitrarily shaped conductor surfaces. We derived an approximate analytical formulation via a perturbative approach consisting in adopting a deformed ellipsoidal coordinate system. We analytically showed the intuitive result that the charge density takes divergent values on the border of two-dimensional conductors providing the closed analytical expression of the quantities of interest. The approach is applicable when, developing the border equation in Fourier series, the first nonconstant border Fourier coefficient is small compared to the constant term, regardless of the geometry of the edge. To the zero-order, we showed that the problem is solved by assigning to the constant A the value π −2 for all the borders. We compared the analytical results with a numerical analysis based on the Finite Element Methods (FEM). For this purpose we analyzed several polygonal contours, and the agreement between the two approaches, already at the zero perturbative order, was found remarkable. This allows us to perform an analytical treatment of the density charge distribution on an arbitrary shaped conducting surface. The presented approach can overcome the difficulty of dealing with divergent quantities when a series approach is used.