A hybridizable discontinuous Galerkin method for simulation of electrostatic problems with floating potential conductors

In an electrostatic simulation, an equipotential condition with an undefined/floating potential value has to be enforced on the surface of an isolated conductor. If this conductor is charged, a nonzero charge condition is also required. While implementation of these conditions using a traditional finite element method (FEM) is not straightforward, they can be easily discretized and incorporated within a discontinuous Galerkin (DG) method. However, DG discretization results in a larger number of unknowns as compared to FEM. In this work, a hybridizable DG (HDG) method is proposed to alleviate this problem. Floating potential boundary conditions, possibly with different charge values, are introduced on surfaces of each isolated conductor and are weakly enforced in the global problem of HDG. The unknowns of the global HDG problem are those only associated with the nodes on the mesh skeleton and their number is much smaller than the total number of unknowns required by DG. Numerical examples show that the proposed method is as accurate as DG while it improves the computational efficiency significantly.


INTRODUCTION
Isolated conductors exist in a wide range of electrical and electronic systems, such as electrode cores of high-voltage inductors 1 , metallic separators of IEC surge arresters 2 , defects in ultra-high-voltage gas-insulated switchgear 3 , passive electrodes of earthing systems 4 , conductors of floating-gate transistors 5 , and, more recently, metallic nanostructures extensively used in optoelectronic devices 6 .In electrostatic simulations of these systems, these conductors result in equipotential surfaces with unfixed (i.e., floating) electric potential values (which depend on the simulation parameters and the geometries of the structures involved) and are referred to as floating potential conductors (FPCs).
Even though execution of these electrostatic simulations by using a finite element method (FEM), which solves the Poisson equation, has become a common practice, accurate and efficient incorporation of FPC models within a FEM framework is still not a trivial task.In recent years, various techniques have been introduced to address this challenge 7,8,9,10,11,12,13,14 .The most commonly used methods among these are the virtual permittivity method (VPM) 7 , the matrix reduction method (MRM) 10 , and the charge simulation method (CSM) 8,9,13 .Each of methods has pros and cons regarding the accuracy and efficiency of the solution, the ability to account for charges on FPCs, and the ease of implementation 11,12,15 .VPM uses a dielectric material with a very high "virtual" permittivity to approximate the conductor.This method is very straightforward to implement since it does not require any modifications to be done on an existing FEM code.But its accuracy depends on the value of the virtual permittivity.Accurate representation of a conductor requires a very high virtual permittivity value but this, in return, makes the FEM matrix ill-conditioned and leads to a less accurate solution.MRM produces more accurate results but it requires rather significant modifications to be done on the original FEM code 10,11,12 .In addition, when an FPC is charged, a nonzero charge condition must be imposed on its surface.Both VPM and MRM can not account for this nonzero charge condition 10,11,12 .
CSM can account for charge conditions since it enforces a specific charge distribution on an FPC but this requires a priori knowledge of simulation results or multiple iterative simulations 10,11,12,13 .Several boundary element methods (BEMs) have also been developed for modeling FPC in electrostatic simulations 2,16 .In 2 a total electric charge condition is applied to determine the potential of uncharged FPCs.BEM is often preferred over FEM for unbounded problems with homogeneous or piece-wise homogeneous materials.In 16 , the Poincare-Steklov operator is used to enforce constraints corresponding to the floating potential.
Recently, it has been shown that FPCs can easily be accounted for using the discontinuous Galerkin (DG) method 15 .By weakly imposing a so-called floating potential boundary condition (FPBC) through the numerical flux, DG can accurately model FPCs with non-zero charge conditions.The implementation of this approach in existing DG codes is rather straightforward.In addition, by enforcing FBPC on the surfaces of FPCs, the requirement to discretize their volumes is removed, which reduces the total number of unknowns.However, even with this reduction, the number of unknowns required by DG is still larger than that of the traditional FEM, which might lead to a considerable increase in computational cost depending on the problem being analyzed.
In recent years, the hybridizable DG (HDG) 17 method has drawn a lot of attention.HDG addresses the fundamental weakness of DG, i.e., reduces the total number of knowns by applying a static condensation technique within the DG framework 18,19 .In HDG, a hybrid variable is introduced on the mesh skeleton and mesh elements exchange information only through this variable.
This approach allows for locating globally coupled degrees of freedom only on the mesh skeleton and results in a global matrix system (in the unknown hybrid variable) with a dimension much smaller than that of the matrix system generated by DG.The local unknowns are then recovered using the hybrid variable that is obtained by solving this smaller global system.Furthermore, HDG can achieve superconvergence by applying a local post-processing technique where the solution represented using polynomials of order leads to a convergence of order + 2 in accuracy.HDG has been applied to various problems 18,19 , such as fluid dynamics 20,21 , and acoustic, elastic, and electromagnetic wave propagations 22,23,24 .Comparative studies 25,26 between HDG and FEM show that HDG outperforms FEM in efficiency for some cases, while it retains other advantages of DG 25,26,18 .
In this work, an HDG-based framework is proposed to implement FPBC in electrostatic simulations.The local problem is formulated with a Dirichlet boundary condition and the electric potential is chosen as the hybrid variable in the global problem.
The floating potential values on each FPC are also left as unknowns in the global problem.The dimension of the resulting global system is equal to the number of nodes on the mesh skeleton (where the hybrid variable is defined) plus the number of FPCs, which in total is significantly smaller than the dimension of the matrix system that would be solved by DG.Other advantages of using the FPBC, such as the ease of implementation, the high solution accuracy, and the ability to account for non-zero charge conditions are inherited by this HDG-based framework.Table .1 compares the properties of this proposed method to other methods briefly described above.Note that the low efficiency of DG is because of the larger number of unknowns, while for CSM it is due to multiple-simulation requirement 12 .
The rest of the paper is organized as follows.Section II starts with the mathematical model of the electrostatic problem in the presence of FPCs, then it presents the HDG formulation for this problem, including the hybridized strong form, the weak form, and the discretized matrix system.Section III presents several numerical examples and Section IV provides a summary.

Mathematical Model
Consider the electrostatic problem described in Figure 1.
In ( 1)-( 4), ( ) is the electric potential distribution to be solved for, ( ) is the permittivity, ( ) is the charge density, ( ) and ( ) are the coefficients associated with the Dirichlet and Neumann boundary conditions, respectively, and ̂ ( ) denotes the outward normal vector of the corresponding surface.Equation ( 4) represents the physical conditions on FPCs.On each FPC, the equipotential value is an unknown, and the total charge is assumed known.The charge condition in (4), i.e., the total electric flux is equal to the total charge, provides the constraint that makes unique.

The strong form
To develop the HDG method, (1)-( 4) are expressed as a first order partial differential equation system by using the electric field Hereinafter, the explicit dependency on is omitted for brevity of the notation.Following the hybridization approach developed in 17 , the local problem on element is defined as where , , and ̂ are the local variables defined on element , ( ) = , is the permittivity in element , which is assumed constant in the element.̂ and ̂ are hybrid variables, which satisfy a global problem (to be described below).
In particular, ̂ is the value of the floating potential on the boundary Γ .Without loss of generality, multiple ̂ variables can be defined for multiple FPCs independently.Equations ( 10)-( 15) describe a local BVP on each element.Once ̂ and ̂ are known, they can be used as Dirichlet boundary values to solve this local BVP.Note that this also means the local variables and can be expressed as functions of the hybrid variables ̂ and ̂ .
The hybrid variables are required to satisfy the global problem defined using the transmission condition 17,18 ̂ ⋅ ( ) = 0, on Ω ∩ Γ (16)   and the charge condition 27 respectively, where ⊙ = ⊙ + + ⊙ − defines the jump at the interelement boundaries and is the total number of faces on the FPC.Since and are functions of the hybrid variables, ( 16)-( 17) can be cast into a global matrix system with unknowns ̂ and ̂ .
The hybridized system described above is amenable to the static condensation of continuous finite element methods (CFEM) and the hybridization of mixed finite element methods (MFEM) 17,18 .Different from CFEM and MFEM, HDG solves the local BVP by DG and enforces the transmission boundary condition weakly using the numerical flux 17,18 .This gives rise to the generalization of the FPBC from DG, where the FPBC is weakly imposed using the numerical flux 15 , to HDG.

The weak form
Let ℙ denote the space of polynomial functions of degree at most ( ≥ 1), then the following discrete finite element spaces can be introduced Let (⋅, ⋅) Ω denote the 2 inner product in the domain Ω ( , and ⟨⋅, ⋅⟩ Γ denote the 2 inner product on the face Γ Following the classical DG approach 28,17,15,27,29 , the weak form of the local problem (10)-( 15) is defined as where ∈ , ∈ , and ̂ ∈ are the approximate solutions sought for (for the sake of simplicity, the same notations are used for the variables in the weak form and the strong form), and the weak form of the global problem ( 16)-( 17) is defined as where the numerical fluxes * and ( ) * are chosen as The numerical fluxes in (22) follow from the local DG (LDG) method 30 and the resulting HDG method is also called LDG-H method 17,18 .The stabilization parameter is of order 1∕ℎ 17,18 , where ℎ is the element edge length.

The discrete system
Choosing the shape functions as Lagrange polynomials 28 , the nodal interpolation of and in each element and that of ̂ on each face are where , = 1, … , , and ̂ , = 1, … , , are Lagrange polynomials, and are the number of nodes per element and per face, respectively, and , , and ̂ are the nodal values.
Applying Galerkin testing to ( 23)- (26) gives the following matrix system where the local unknown vectors are defined as the global unknown vector is defined as and ̂ is the local copy of the global unknowns, i.e., each local face of element is mapped from a face of Γ Here, = 1, ..., , = 1, ..., .Fig. 2 illustrates the mapping between the nodes of local elements (blue dots) and the nodes of the skeleton (red circles).
Note that has dimensions × , where is the size of the operand (input vector) and is the size of the output vector .
Solving [ , ] in terms of ̂ and ̂ from ( 27) yields where Inserting ( 32) into ( 28) yields a global system involving only the global unknowns ̂ and where The dimension of the global system (34) is ∼ ( + 1), which is much smaller than that of the DG method, which is ∼ , as described in 15 .Once ̂ and ̂ are solved from (34), they can be used to solve [ , ] in the local system (32).
Since the local problems of different elements are independent from each other, they can be solved in parallel.As the dimension of ( 32) is only ∼ , the computational cost of this step is relatively low and can be ignored, especially in large scale problems 18 .
Figure 3 (c) illustrates the nodes where and ̂ are defined.Here, the case of = 5 is illustrated.The degrees of freedom of ̂ only correspond to the nodes on the wireframe while those of correspond to all nodes.Same as DG, the nodes of are doubly defined on the wireframe.Table .2 compares the dimensions of the global matrix systems of DG 15 and HDG.One can see that the difference in matrix dimensions becomes more significant with increasing order of the polynomial basis functions.This is easy to understand from Figure 3 (c).As the polynomial order is increased, the ratio of the number of nodes in the interior of each element to that on the element surface increases.

Plasmonic-enhanced Photoconductive Antenna
Next, a plasmonic-enhanced photoconductive antenna (PCA) is considered.The device geometry is shown in Fig. 4 (a).The semiconductor layer (blue) is made of GaAs, which is a photoconductive material that can absorb optical electromagnetic (EM) wave energy and generate terahertz (THz) signals 29 .The metallic nanostructures (yellow) are designed to enhance the local EM fields and hence increase the optical-to-THz efficiency 6 .Throughout the operation of this device, a bias voltage is applied on the two electrodes (red), generating a static electric field.Thus, to model this device, one needs to solve the electrostatic problem under the bias voltage, aside from the transient EM response 27,15,29 .
The semiconductor layer has relative permittivity 10.9 and the surrounding area (gray) is air.The computation domain is truncated at the outmost boundaries with homogeneous Neumann boundary condition.Dirichlet boundary conditions are applied on the surfaces of the electrodes, with ( ) = 0 on the left one (cathode) and ( ) = 10 V on the right one (anode).Since the metallic nanostructures are isolated conductors, they act as FPCs and an independent FPBC is applied on each block of the nanostructures.and DG 15 (using the same polynomial order = 4).The maximal difference between the solutions obtained from the two solvers is 1.1×10 −4 V. We also note that no volumetric meshes are used in the FPCs and the electrodes since they are treated as boundary conditions in our method.In practical device simulations, this treatment can save considerable computational resources since finer meshes are usually required near the nanostructures 27,29 .Table .3 shows the dimensions of the global matrix system in DG and HDG.Again, we can see that as increasing the order of the polynomial basis functions, the saving of the number of unknowns from HDG becomes more significant.

Surge Arrester
Next, the proposed method is used to compute the electric potential on an IEC surge arrester 32 .The model is shown in Figure 5 (a).The arrester consists of three segments of metal-oxide varistor (MOV) column.Each segment is surrounded by a porcelain layer.The pedestal and the surrounding cylinder are grounded ( = 0) and a high voltage ( = 100kV) is applied on the lead and the grading ring 32 .The MOV columns are separated by two isolated metal flanges, which are considered as FPCs and are modeled with two independent FPBCs in the proposed HDG method.
The diameters of the MOV, the inner wall of the porcelain layer, the outer wall of the porcelain layer, the flanges, the pedestal,  17kV.These results agree with the data reported in 32,16,2 .The maximal difference between the solutions of HDG and DG over the whole domain is 1.1 × 10 −5 kV.The dimensions of the global linear systems of HDG and DG are 791, 282 and 1, 127, 672, respectively.

CONCLUSIONS
A HDG scheme for modeling FPCs in electrostatic problems is developed.The local problem is formulated as a Dirichlet BVP, the global problem is formulated in the unknown electric potential and the unknown floating potential values of each FPC.
The proposed HDG scheme retains the advantages of the DG scheme previously proposed for FPC modeling in electrostatic simulations, i.e., FPCs that can account for non-zero charge conditions, accurate solution, and ease of implementation in an existing code.Meanwhile, it significantly reduces the number of degrees of freedom as compared to DG, leading to a reduced computational cost.

FIGURE 1
FIGURE 1 Schematic description of an electrostatic problem involving multiple isolated conductors.

FIGURE 2
FIGURE 2 Illustration of the local unknowns (blue dots) and the global unknowns (red circles) in HDG with = 3.
(a).A thin metal tube is inserted into a coaxial capacitor.The voltages applied on the inner and outer boundaries of the capacitor are (| | = 0 ) = 0 and (| | = 1 ) = 1 , respectively.The metal tube is modeled as an FPC and the FPBC is applied on | | = 2 and | | = 3 .The total charge on the FPC is .The analytical solution of the electric potential is

FIGURE 3
FIGURE 3 (a) Schematic description of the coaxial capacitor model.(b) Comparison of numerical solution obtained by the proposed HDG with analytic solution on the line ( , = 0) for different values of .(c) Illustration of the nodes where and ̂ are defined.

Fig. 4 (FIGURE 4
Fig. 4 (b) shows and ̂ on the plane ( , , = 3 m) solved from the proposed HDG method with = 4.As expected, ( ) is constant on the surface of each FPC block.Meanwhile, the potential values on different FPC blocks are different because theFPCs are isolated from each other.One consequence of the inhomogeneous electric potential distribution is that strong local static electric fields are generated near the FPC, which greatly influences the carrier mobilities in the semiconductor layer27 and hence influence the device performance31 .
and the lead are 60 mm, 140 mm, 200 mm, 230 mm, 280 mm, and 40 mm, respectively.The major and minor diameters of the grading ring are 1130 mm and 70 mm, respectively.The heights of the pedestal, each metal flange, each MOV segment, and the lead are 2000 mm, 120 mm, 960 mm, and 3400 mm, respectively.The diameter and height of the surrounding cylinder are 8000 mm and 9000 mm, respectively, which are determined by the minimum phase-to-earth clearance32 .The relative permittivity of the MOV column and the porcelain layer are = 800 and = 5, respectively.

TABLE 1
Comparision of FEM-based methods used for modeling FPCs in electrostatic simulations.

TABLE 2
Dimension of the global matrices in DG and HDG for the coaxial capacitor example.

TABLE 3
Dimension of the global matrices in DG and HDG for the plasmonic-enhanced PCA example.