Plate capacitor problem as a benchmark case for verifying the finite element implementation

In this work, parallel plate capacitors are numerically simulated by solving weak forms within the framework of the finite element method. Two different domains are studied. We study the infinite parallel plate capacitor problem and verify the implementation by deriving analytical solutions with a single layer and multiple layers between two plates. Furthermore, we study the finite parallel plate capacitor problem and verify it by Love’s potential equation and Xiang’s capacitance equation. Moreover, the fringing effect is considered and extended to problems with multiple dielectric layers, such a solution is not possible by means of the existing analytical solutions. Besides, we realize the possibility of choosing different boundary conditions (electric potential boundary conditions and charge density boundary conditions) by changing the weak form. Finally, a transient solution that includes dielectric loss and calculates the quality factor of a capacitor is presented, which may be used in capacitor design. Convergence and consistency of results are demonstrated by comparing the results between analytical and numerical solutions and also the results from different boundary conditions.

A parallel plate capacitor problem with given dimensionless coordinates, ρ, ζ, ζ and a spacing parameter, κ. Drawing is taken from [2,Fig. 1] equation. This so-called Love's equation is solved in cylindrical coordinates in [4]. The electric potential, V , in Love's theorem is given by where both square roots have positive real part and f (t) is the unique solution of the integral equation These solutions have been presented in [3]. The parameters in Eqs. (1) and (2) refer to Fig. 1. Equation (2) is a Fredholm integral equation of the second kind and various numerical techniques have been applied to solve this integral equation [5][6][7][8]. In [9] a different analytical solution for parallel plate capacitor problem is proposed, which is criticized in [10] since it fails to satisfy all the requirements of the problem. The electric field outside the parallel plate capacitor in [11] may be captured by an explicit analytical solution as well. By reducing the system, the distance between the two plates has been simplified to vanish, such that the solution represents the electric field distribution outside the plates. A generalized Love's equation is motivated in [12] for a capacitor composed of two discs of different radii.
Capacitance is one of the basic concepts in electrodynamics measuring the performance of a capacitor by a ratio between charge and potential difference of separated plates. For example, a capacitor constructed by two parallel plates, both of area A, separated by a distance d, an "ideal" capacitance is given by where ε 0 is vacuum permittivity. Under the condition A d 2 , this simplified solution holds accurately. The case in Eq. (3) holds true whenever the area between plates is under vacuum or filled with air. When a parallel plate capacitor has a dielectric with a relative permittivity, ε r , between the plates, then the capacitance reads These ideal capacitance equations assume that the electric field is uniform and also perpendicular to the capacitor electrodes. These simple forms fail to account for fringing effect (also called edge effect). The fringing effect represents the non-uniform electric field around the edge of a parallel plate capacitor (see  2). When studying a finite parallel plate capacitor problem, edge effects need to be considered in order to improve the accuracy of the capacitance results. Kirchhoff studied a circular parallel plate capacitor problem in 1877, and gave an equation for the capacitance with fringing effect [13]. Kirchhoff's equation has been proved after 85 years in [14]. Technically, the determination of the capacitance of a plate capacitor is directly related to Love's equation. In order to obtain the capacitance, different numerical and analytical approaches have been developed [15][16][17]. One often used technique is the small separation assumption, which is based on the simplification that two plates are separated by a very small distance regarding their lengths. A more general approach in [18] characterizes the capacitance for larger separation per length ratios. In [19] Xiang proposed a capacitance equation as shown in Eq. (5) that takes into account fringing effect for two-dimensional parallel-plate capacitor problem, which is used to verify a finite element method based numerical simulation in [20]. The λ coefficient in Eq. (5) is related to L/d. Because of the two-dimensional problem, L represents the plate length instead of the area, i.e. the length along the direction perpendicular to the two-dimensional plane, which is regarded as the unit length such that the area is L × 1. In this paper we use Xiang's equation to verify our finite element results of two-dimensional parallel-plate capacitor problem and we explain the following result in detail, The problem of a dielectric between two plates is challenging. In [21] the dynamics of the terraces on the silicon surface due to the combined effect of strain and current are investigated. In [22,23] the necessary material parameters appearing in the thermo-electric-viscoelastic constitutive model are determined experimentally and thus the relevant numerical model is developed. In some cases, for example by replacing the solid dielectric by an insulating fluid flow, the dielectric even varies between the two plates. The design of multilayer dielectrics has found application in the study of dielectric elastomer actuators [24][25][26][27]. Yet the literature lacks for studies on the problem of single/multiple dielectrics between two plates, especially in multiple dielectrics, albeit their use in electronics is quite common. For solving electric field, we use governing equations with the electric potential as a given boundary condition. As a consequence, two plates of the capacitor store equal amounts of opposite charges that create the electric potential difference and an electric field. Theoretically, the charge density may be used as a given boundary condition for such a study as well. We address the question of differences in choosing electric potential versus charge density as the boundary condition. Especially in functionalized surfaces, where the surface charge is manipulated by a coating, it is of utmost importance to comprehend the physics related to the charge and electric potential, where they are inherently coupled by Maxwell equations [28]. In reality, electrical components often work under a periodic electric field environment, which results in energy dissipation in a capacitor. At present, there are almost no literature discussing the numerical analysis of this problem.
In this work, we develop the weak form of the governing equations and solve it by the finite element method and verify the implementation by deriving an analytical solution with a single/multiple layers. The finite element model is developed by open-source packages, the FEniCS project. 1 All codes are made public on Github 2 to be used under the GNU Public license 3 for promoting an efficient scientific exchange as well as further studies. We demonstrate a verification of the numerical implementation across the length-scales, since the Maxwell equations are valid for all length-scales in continuum mechanics. By choosing boundary conditions differently, we expect building models for selectively studying various combinations of materials for dielectric and plates leading to a better understanding of the underlying mechanisms. Additionally, we calculate a capacitor's quality factor and give a transient solution that takes into account dielectric loss, which is modeled by an irreversible part of charge potential.

Governing equations
We begin with classical electrodynamics theory, we refer to [29][30][31][32] for a more detailed introduction. The Maxwell equations are introduced in two groups, from Faraday's law, we obtain where B is the magnetic flux (area) density in T (tesla), E denotes electric field in V/mm (volt per meter). We use standard continuum mechanics notation with ∇ for grad operator meaning a space derivative and ∂()/∂t as the rate. A solution of these equations is established by the following ansatz functions: with the electric potential, ϕ in V, and magnetic potential, A in T mm. This solution introduces gauge freedom, in other words, we may select ∇ · A and ∂ϕ/∂t freely without changing the electromagnetic fields, E, B. From the balance of electric charge, with the charge (volume) density, q in C/mm 3 (coulomb per volume), we obtain the following Maxwell equations: where the charge potential, D in C/mm 2 , current potential, H in mA/mm (ampere per meter), and electric current (area) density, J in mA/mm 2 , need to be defined by constitutive equations. The form of Maxwell equations depend on the chosen unit system and we use the Maxwell-Lorentz aether relations: where ε 0 = 8.85 × 10 −12 mA s/(V mm) and μ 0 = 4π × 10 −13 V s/(mA mm) are universal constants. The parallel plate capacitor problem is an electrostatics problem such that magnetic potential vanishes, A = 0. The total charge is composed of free and bound charges, where the free charge potential, D, and bound charge potential, P, are introduced, as follows: where P is also known as electric polarization. For a linear material, we use where dielectric permittivity, ε r , is a material parameter. Therefore, from Eq. (8), we obtain In Eq. (12), q f is the free charge per unit volume, which is different from surface charges. The electric potential is a (smooth) continuous function, but the charge potential varies for different materials (mediums) with differing permittivity, ε r . The top and bottom boundaries of the dielectric are such singular interfaces, where the charge potential experiences a jump across these interfaces, see the top and bottom boundary of blue area in Fig. 3. We assume that there are no free charges on the interfaces, q f I = 0 and obtain where n is the unit normal vector of interface and · indicates a jump between values at either sides of the interface, D ± . Equivalently, by using Eq. (11), it reads n · ε 0 ε r ∇ϕ = 0, Within the dielectric, no free charges exist, such that motion of free charges and thus rate of free charge potential vanish. In other words, the charge potential is instantaneous and has to satisfy Into the latter we insert Eq. (11) and obtain the well-known Laplace's equation for the electric potential,

Analytical solution
We consider an infinitely long parallel plate capacitor problem in two-dimensional continuum as illustrated in Fig. 3. The red lines are the top and bottom plates made of a conductor (metal). The origin of the coordinate system is located at the center of the capacitor. The area with the blue shadow is a dielectric material and the white dashed line is the centerline of the dielectric. The geometry is defined by the distance between the two plates, H , the thickness of the dielectric, h, and the offset of the dielectric center relative to the x-axis, c. The relative permittivities of two different media are ε r1 and ε r2 . Aligned with the infinitely long capacitor along x, we assume that the problem is reduced from a two-dimensional problem to a one-dimensional problem by letting the electric potential, ϕ, vary only along y-direction. This simplification results in the following boundary value problem: The Laplace problem in Eq. (16) reads By integrating twice, we obtain with the gradient C a andD a are coefficients. At the top interface, y = h 2 + c with n = e y , we have as a limit value ε + r = ε r2 and ε − r = ε r1 . With the solution in Eq. (20), we obtain (∇ϕ) + =Ĉ 2 and (∇ϕ) − =Ĉ 1 . Inserting these values into Eq. (14), the transition condition at the top interface reads ε r2Ĉ2 = ε r1Ĉ1 . Analogously, for the bottom interface, y = − h 2 + c, we obtain ε r2Ĉ3 = ε r1Ĉ1 . We summarize the transition conditions, In order to fulfill the Laplace's equation (16) the electric potential, ϕ, has to be continuous in the entire domain. Therefore, on the interfaces, we obtain the continuity conditions, By using them in Eq. (19), the continuity conditions read As boundary conditions, we define the electric potential as given on both plates, By using the transition conditions in Eq. (21), continuity conditions in Eq. (22), and boundary conditions Eq. (23), we obtain for the electric potential, We continue in order to reach the goal of getting the capacitance of this system. For the problem depicted in Fig. 3, where between plates different dielectric materials (air and dielectric material) are used, we may macroscopically divide the entire parallel plate capacitor into a parallel/series problem of several capacitors composed of a single dielectric, as shown in Fig. 4. In the case of a parallel connection in Fig. 4a, the capacitance is calculated by In the case of in series connection as in Fig. 4b, the capacitance reads We propose to use these simplified calculations above, since the Laplace problem is a linear differential equation and the superposition principle holds. These results are used to verify the numerical solutions in the following.

Computational implementation
The finite element method (FEM) is used by following implementations in [33,34] and refer to [35,Sect. 3] for engineering examples in electromagnetism by using FEniCS packages. Herein, we implement two FE models: the finite parallel plate capacitor in Fig. 5a and the infinite parallel plate capacitor in Fig. 5b; both of them as two-dimensional problems with upper and lower plates treated as interfaces. The computational domain, , is composed of 1 and 2 representing two different mediums, = 1 ∪ 2 . The closure of the domain, ∂ , is the boundary to be defined by additional conditions. For the finite plate problem, a sufficiently large domain is needed to ensure that the FEM corresponds to the analytical solution. In the case of the infinite plate problem, periodic boundary conditions are used on left and right sides reducing the size of the computational domain significantly. In other words, the solution space is constructed with a restriction that the solution on left and right ends are identical. Thus, the mesh needs to be utilized in such a way that the vertical coordinates of nodes on the left are matching the vertical coordinates on the right. Periodic boundary conditions are used to model infinity only along the x-direction in Fig. 5b. Boundaries are denoted by 1 , 2 . . . 7 . By using a variational formulation, we begin with the governing equation (15), multiply by an arbitrary test function, δϕ, and obtain an integral form within the domain composed of finite sized domains E called elements, The electric potential, ϕ, is the unknown variable and the constitutive equation (30) involves first gradient of ϕ already, The integral form incorporates second-order partial derivatives of ϕ, whereas no derivatives of δϕ. According to the Galerkin procedure, we may use the same space for ϕ and δϕ such that we choose to weaken the continuity condition by integrating by parts and applying Gauß's theorem, We emphasize the used notation; • Between elements, I represents the interface used herein modeling the upper and lower plates. The sum over elements introduces the jump brackets since neighboring elements have oppositely directed surface normals, n, on the shared surface. • The sum over elements let the terms on shared boundaries vanish since D is continuous except on the interface. • With a sum over elements E ∂ E \I = ∂ represents Lipschitz boundaries, which is sufficiently regular and possesses only outer boundaries. • Within the domain, we obtain after assemble, E E = .
In the Sect. 2, we assumed that the charge on the interface is 0, where the interface referred to the interface between different media. But in this section, I refers specifically to two plates, where we have no zero charge assumption. Therefore, how the second term in the Eq. (31) is handled depends on different boundary conditions, which is explained in detail below. The standard FEM elements are employed [36], they are of order q and span P q on E . Shape functions are piecewise or elementwise continuous over the whole domain with a compact support on the element such that nodal values have influence only locally. This representation is adequate for approximation in Hilbertian Sobolev space H 1 . Piecewise continuity ensures a monotonic convergence. Linear Lagrange elements are used for ϕ, which means we use interpolation polynomial of first order as shape function for each element in scalar space, In addition, because of the existence of dielectric, D is piecewise continuous in vector space. We choose discontinuous Lagrange elements for D, which implies L 2 elements with P 0 such that they generate a constant within the element with a jump across the boundaries. We emphasize of not using any special elements [37] designed for electromagnetism herein. Therefore, an extension in multiphysics as in [38,39] is straightforward that is a benefit of the presented formulation. In Eq. (31), the integral over the interface, I , is related to the charge density on the interface. For implementation of this term, we introduce two distinct cases: A: We use the electric potential, ϕ, as given on the two plates (interface). This so-called Dirichlet boundary condition let the boundary term vanish. since δϕ = 0 is set on Dirichlet boundaries. B: Instead, we may apply a charge density as a Neumann boundary condition and observe the significance of this term.
Next, we will clarify these two boundary conditions (case A and case B) separately.

Case A: Electric potential boundary condition
As we use Dirichlet boundary condition on the interface, its solution is given such that the nodes are taken out of the system, in other words, the numerical implementation circumvents from solving the values on the corresponding degrees of freedom. The weak form (31) simplifies to ∂ δϕD · nds − D · ∇δϕdv = 0.
For the model as depicted in Fig. 5a, on 1 we set ϕ 1 = 0 to simulate zero potential at infinity-this boundary, 1 , denotes the boundary of the whole computational domain. Hence, we need a sufficiently large domain in order to incorporate this far-away boundary condition correctly. 2 is the boundary between two different media, 1

Case B: Charge density boundary condition
The charge stored on the plate is the source of the electric potential. In order to introduce the charge density on the two plates as the boundary condition, we need to change the second part in Eq. (31) and introduce the following relationship: where q f is the charge density. We fail to obtain q f directly from Eq. (34) in FEM, but we may easily calculate the integral over the interface in order to obtain the total charge: Q is the total charge on one plate and another plate carries the same amount of negative charge. For having cases A and B identical, the total charge, Q, is calculated from the results acquired by solving Eq. (33) in Case A.
Substituting Eq. (34) into the second part of Eq. (31), we get the weak form for Case B, For the charge density boundary condition, we simulate Fig. 5b. Because charge distribution on an infinite parallel plate capacitor is uniform, the charge density is simply found by q f = Q/L, where L = I dS is the length of 5 and 6 . The charge density boundary conditions read

Results
Capacitance is an important physical quantity for parallel plate capacitors and it is also one of the results in focus. We already have elementary relations on Eq. (27) and Eq. (28) to get reference values of capacitance. In FEM we need to calculate capacitance from the electric field results that we have solved. First, we introduce the electric field energy. It is the energy contained in an electric field. The energy density or energy per unit volume of an electric field becomes The total energy U E stored in electric field with the domain, , is In addition, the total electric energy stored in a capacitor is given by where C is the capacitance and V is the electric potential difference between the two plates. Using Eqs. (39) and (40) we get the final equation: With the electric field distribution results solved by the FEM, we easily calculate Eq. (41) for obtaining the capacitance, C. Taking model (a) as an example, we select the volume between the two plates as the calculation area. As shown in the Fig. 6, Eq. (41) becomes  Table 1 Parameters setting for the analytical solution, L-length of the plates. In the FE model (b) (Fig. 5b)

the domain size along y-direction A = 50
Parameters Values

Comparison between analytical and numerical solutions
Infinite parallel plate capacitor model in Fig. 5b has been studied with the parameters as compiled in Table 1 referring to Fig. 3. We compare the capacitance computed by the numerical solution using different mesh sizes to the analytical solution. We use degrees of freedom (Dofs) in order to represent the size of mesh. As known as the hconvergence study, we increase degrees of freedom by preserving the element quality and the result in Fig. 7 indicates that the numerical solution is monotonously converging to the reference solution that is herein the analytical solution from Eq. (28). The difference is caused by Eq. (28) since we use different dielectrics and  25) and (26). Numerical results by solving Eq. (33) series connections for the estimation. The comparison over one value may be misleading, therefore, we analyze the solution in the whole domain. We compare the electric potential and electric field results between two plates at x = 0 along y-axis. In Fig. 8 we observe that the numerical solution results match the analytical solution in the whole domain.
For the finite parallel plate capacitor problem, model (a), we compare the results herein with a numerical result that is presented in [6] of Love's potential equation (1). We study two different capacitors regardless of the dielectric with L/ h = 2 and L/ h = 20 (L-length of the plates, H -distance between the two plates) according to [6]. The comparison of the equipotential lines results is shown in Fig. 9. As shown in Fig. 9, we specify the same level of equipotential lines as in [6]. The distributions of equipotential lines obtained by the two different methods are almost identical.
For the capacitance, we compare FEM result with Xiang's equation: and where L is the length of the plate, h is the distance between the two plates. Equation (44) implies that λ < 1/2 because 0 < tanh (π L/2 h) < 1. The power series on the right side of Eq. (43) converges very quickly and we have chosen to include four terms, as follows: The comparison of our FEM capacitance results and capacitance from Eq. (45) is shown in Fig. 10. On the x-axis, we use L/ h representing the ratio of the plate length to the distance between the two plates. The larger L/ h means that the parallel plate capacitor becomes flatter and we observe that the two results agree better. We conclude to have verified the accuracy of the finite element model proposed in this paper. Both the potential and capacitance results show good agreement with the analytical solutions for both infinite and finite parallel plate capacitor problems. It should be pointed out that for the finite plate problem neither Love's equation nor Xiang's equation considers the presence of a single layer or multiple layers between two plates. In other words, the above equation only applies when the material between the two plates is homogeneous and its relative permittivity is ε r = 1. We may calculate the capacitance of a capacitor with a single layer or multiple layers according to Eqs. (27) and (28). As mentioned before, these two simple equations do not account for edge effects, so the results are less accurate for the finite parallel plate capacitor problem. For this kind of problem, FEM shows good adaptability and flexibility, since the proposed implementation not only allow to study the presence of various dielectric layers between parallel plate capacitors, but also have the potential to study the dielectric layers with irregular shapes. Next we study the capacitance of parallel plate capacitors with multilayer dielectrics and use the analytical results of Eqs. (27) and (28) as reference.  For multilayer dielectrics, we study two cases, where the difference between them is the diverse distribution of dielectrics in the capacitor as shown in the Fig. 11. For these cases we set the size of domain A = 20 m referring to Fig. 5a. This size is large enough compared to the size of the parallel plate capacitor in the two cases. And the capacitance results are shown in Table 2. Numerical and analytical results are quantitatively the same. When compared with the analytical results, the numerical capacitance results obtained by FEM used in this paper are larger in Table 2. We understand that the FEM considers fringing effects correctly where it is expected to obtain larger results because of the electric field concentration around the edges. The difference between the FEM and ideal capacitance results of Case 2 is slightly larger than that of Case 1, because the H/L of Case 2 (H/L = 0.5) is larger than Case 1 (H/L = 1/6), and the fringing effects increase as H/L increases. We visualize the electric field for presenting the fringing effects in Fig. 12, where electric field magnitude distribution indicates the concentration at the ends of the plates.

Comparison of two different boundary conditions' implementations in the multilayer approach
We use a similar configuration as before, but this time a multilayer setting as the geometry is shown in Fig. 13. We choose a A d such that the simulation corresponds to an infinitely large plate along the y-direction. Parameters are compiled in Table 3.
Initially, we solve Eq. (33) under electric potential boundary conditions in order to obtain the electric potential and electric field results. Then, we calculate the total charge, Q, on one plate by means of Eqs. (34), (35). We use this total charge result as the value of Q in Eq. (37). In spite of, we use different boundary conditions, they are theoretically equivalent and the two results should be the same under ideal conditions.
The results are shown in Table 4 and Fig. 14. We observe that the difference between the results of two boundary conditions is around 1% of relative error. We stress that the charge density BD includes a numerical rounding error because of calculating Q and 1% error is reasonable to claim that the error is caused solely by this fact.

A transient solution by including dielectric loss
The computational implementation is generic to be used in three-dimensional continuum as well as ready for parallelization by using mpirun. Therefore, industrial applications may be simulated as well. We construct the C2220 capacitor in Fig. 15, where the material of two plates are made of copper and the relative permittivity of copper is ε Cu r = 1 since it is an electric conductor. We choose PTFE (commercial name is Teflon given by DuPont de Nemours, Inc.) as dielectric with relative permittivity ε PTFE r = 2.1. We simulate a transient response by a harmonic boundary condition, where the maximum amplitude distribution may be depicted in Fig. 16.
For time varying electric fields, flux of the electric energy is reversed leading to a necessity of dielectrics also as mechanically supports in separating electrical conductors. The mechanical deformation of dielectrics creates losses in this energy transfer modeled as propagating waves. A possible model relies on permittivity, The purpose is to model the electric field as a harmonic function by using the Euler equation, where E amp is electric field amplitude. Equation (11) is rewritten as by employing the rate, This material model is clearly considering the reversible part as before as well as a rate term that incorporates losses in a cyclic electric loading. We may rewrite by separating reversible and irreversible (dissipative) terms where e 1 = ε r and e 2 = ε r /ω. For the computational implementation, we apply the finite difference method, also known as Euler-backward approach, where time step, t, is set constant for the sake of simplicity. One time step before, the electric potential is calculated, ϕ 0 , which is then used for obtaining E 0 . By using Eq.(51) for D in Eq. (33), we achieve a simulation with choosing the time-varying electric potential boundary condition ϕ 5 = ϕ 0 2 sin(2πt) and ϕ 6 = − ϕ 0 2 sin(2πt) by referring to Fig. 5b and solving for different values of e 2 in one circle, we plot the hysteresis curve in Fig. 17. Indeed, these dielectric losses are unwanted in design. Therefore, a quality factor is used in a capacitor, also known as the Q-factor, representing the efficiency of a given capacitor in terms of energy losses Herein we denote the quality factor by Q c to distinguish it from the total charge notation. Quality factor Q c of a parallel plate capacitor is then studied using Teflon as the dielectric. In the microwave region at a working frequency of ω = 10 GHz, Teflon has an approximate real part of the relative permittivity ε r ≈ 2.03 and a approximate imaginary part ε r ≈ 0.0008 according to [41]. We choose the time-varying electric potential boundary condition ϕ 5 = ϕ 0 2 sin(ωt) and ϕ 6 = − ϕ 0 2 sin(ωt). We define Q c by Eq. (52), by calculating the reversible energy in one cycle as well as the dissipated energy in this cycle We obtain the quality factor of this parallel plate capacitor, C2220 with Teflon as the dielectric at a working frequency of 10 GHz as Q c = 1652. Obviously, dissipation in each cycle means an energy loss and a capacitor storing an energy has a better quality if the dissipation goes to zero meaning Q c → ∞. In industrial applications, datasheet for a component has an indication about this value for a target frequency, for example good quality capacitors reach at 1 MHz over Q c > 10 000 or for 1 GHz over Q c > 1 000. By using a transient analysis, we demonstrate that it is straightforward to calculate this value for different frequencies. We stress that the frequency changes Q c value. On the x-axis, we use E y evaluated in the middle of the computational domain, and on the y-axis, we use D y evaluated at the same position over time

Conclusions and outlook
In this paper we have scrutinized the weak form of an electrostatic problem in the case of surface charges. Surface charge modeling is of importance on interfaces between phases or different materials exchanging electric charges, but a numerical implementation is difficult to find in the literature [42,43]. The governing equations for electrostatic problems are solved by the finite element method (FEM) in order to study the parallel plate capacitor problem. With a choice of a "simple" problem, we have been using an analytical solution to quantify the accuracy of the numerical implementation. The verification has shown the expected monotonous convergence, which is indeed an important feature in the FEM. Even for the case of boundary conditions with a given surface charge, our proposed implementation has this feature such that a posteriori error analysis is possible for more "complex" problems without a closed-form solution.
Additionally, the charge density boundary condition has been verified by comparing with the equivalent potential boundary condition, which is a simple benchmark test for consistency. The two boundary conditions gave essentially the same results for electric field and potential. Therefore, we conclude that the use of continuous and discontinuous Lagrange elements is adequate for this problem.
Finally, we have applied the methodology for an industrial application by using an irreversible part for the charge potential and determined the quality factor of a capacitor and provide a transient solution that accounts for dielectric loss. Accurate computations are of importance for testing the predictive capacity of reduced order systems [44] as well to design sophisticated structures by using novel dielectric designs [45][46][47].