Parallel Finite Element Computation of Time-Varying Ionized Field Around Hybrid AC/DC Lines via Fine-Grained Domain Decomposition

The time-varying hybrid ionized field around HVAC and HVDC transmission lines is a computationally demanding problem due to the coupling of the Poisson’s equation and current continuity equation, as well as the involvement of large matrix in traditional Galerkin finite element method (FEM). In this paper, a fine-grained nodal domain decomposition (NDD) scheme, which enables each sub-domain with only one unknown to be solved independently in a massively parallel fashion, was employed to solve the Poisson’s equation. Meanwhile, an upwind nodal charge conservation (NCC) method is applied to solve the current continuity equation without numerical oscillation at each finite element nodal level. The computation of NDD and NCC can both be vectorized and mapped to massive computational cores and utilize the computing power of graphics processor units (GPUs). The interaction between HVAC and HVDC was solved without the Deutsch’s assumption to guarantee the accuracy, and the wind influence can be considered. With the massively parallel NDD scheme and NCC scheme, both the Poisson’s equation and the current continuity equation were solved at each time-step on GPUs to obtain the transient details of the hybrid ionized field. The performance of the proposed method is tested and compared with commercial software, showing a speedup of 17 times for an 8184-node finite element case with a mean relative error of 0.07%.


I. INTRODUCTION
The advantage of high voltage direct current (HVDC) over high voltage alternating current (HVAC) are presented in many perspectives: higher power transmission capacity, lower net cost for long-distance transmission, no skin effect, lower line loss, and easier system design [10]. Because of these benefits, numerous HVDC-related projects have arisen in recent years [5]. Among these projects, some are retrofitting of existing multi-circuit AC transmission towers with DC lines to increase power transfer because of the high cost of building new HVDC transmission towers, restrictive access to new right-of-way, and long duration transmission planning.
The proximity of AC circuits and DC circuits on the same transmission towers has an appreciable influence on The associate editor coordinating the review of this manuscript and approving it for publication was Tomás F. Pena . the ionized field around them due to corona interactions between two circuits [15]. An efficient numerical computation for the ionized field is therefore desired by engineers for designing the conductor configuration and tower geometries, and assessing environmental impact. The computational efficiency is important for a hybrid ionized field because this is a time-marching problem with a large problem domain which requires repeated FEM computation to handle multiple degrees of freedom (DoF).
Research in the ionized field has been carried on for more than a century. However, computation efficiency wasn't a major focus. In 1914, the first analytical solution for coaxial cylindrical geometry was developed in [23], but it had very limited applications. In 1933, Deutsch laid the cornerstone for numerical analysis by proposing an assumption that the space charges affect the magnitude but not the direction of the corresponding charge-free field. Although Deutsch's assumption simplified the calculation, the validity of this assumption was questioned [8]. Researchers began to explore methods independent of Deutsch's assumption. The first FEM-based iterative method was proposed in 1979 in which the space charge was tuned until convergence using Poisson's equation and current continuity equation [11]. Different methods for solving Poisson's equation and current continuity equation and updating space charge were developed later [14], [21], [28], but the basic iterative algorithm remained, and computation efficiency has not been improved. Moreover, these methods were only applicable to HVDC scenarios because the steady-state result did not fit the hybrid situation.
The ionized field around hybrid transmission lines has been researched for almost four decades. In 1981, Chartier et al. observed that the voltage gradients on the surfaces of AC and DC conductors were time-varying for hybrid transmission lines, which was in opposition to HVDC where the voltage gradient was time-invariant [7]. Then, the significant impact of the AC-DC circuit interaction on the electric field was discovered [15]. In 1992, the finite element method (FEM) was used to iteratively calculate the space charge distribution [3] but the computation was still costly. Deutsch's assumption was used to reduce the cost, however, it reduced the accuracy of the result. A calculation method that decomposed one AC cycle into several DC cases and solved these DC cases using the numerical method for the steady-state ionized field was developed in [27], but the dynamics of the hybrid ionized field were neglected.
However, few materials have been reported to handle the computational efficiency problem encountered in the ionized field computation, especially for the time-varying case. The goal of improving computation efficiency has become feasible with the development of high-performance computing hardware, and a critical step is finding an algorithm that can fully utilize this hardware. In this paper, an algorithm based on FEM for solving the hybrid ionized field is proposed. To obtain the solution for Poisson's equation, a finegrained nodal domain decomposition (NDD) methodology is implemented on graphic processing units (GPUs) to obtain the massive parallelism and hence to improve the computation efficiency. NDD methodology is chosen to provide sufficiently simple sub-domains with only one unknown in each sub-domain so that each sub-domain can be solved independently. The upwind nodal charge conservation (NCC) method is also a part of the algorithm and it is applied to solve current continuity equations. Each finite element node can be projected to a compute unified device architecture (CUDA) core since both NDD and upwind NCC are methodologies based on nodes.
This paper is organized as follows. The assumptions to simplify the problem is given at the beginning of Section II. In the same section, the governing equations and boundary conditions are introduced and explained. Then, in Section III the NDD methodology, upwind NCC method, and application of boundary conditions are described. Next, the hybrid line configuration, the result comparison between the proposed algorithm and COMSOL Multiphysics TM , and the computation speed-up for a hybrid ionized field case are provided in Section IV. Conclusions are drawn in Section V.

II. PROBLEM DESCRIPTION A. ASSUMPTIONS FOR SIMPLIFICATION
To reduce the complexity of corona phenomena around the transmission lines, some reasonable assumptions are employed to build a solvable model [15], [25]: 1) The thickness of the ionization layer around the conductor is small enough to be neglected. 2) Diffusion of positive and negative charge is neglected since it has a very slight effect on the ionized field comparing to the convection.
3) The mobilities of positive ion and negative ion are assumed to be constant. The coefficient of recombination is assumed to be constant. 4) Ions generated by AC conductors due to corona are all restricted in a thin layer around the AC conductor, which is not a part of the analysis domain. 5) The Kaptzov's condition will be guaranteed, which states that the gradient of electric potential on the conductor will not exceed the onset initial gradient when corona happens.

B. GOVERNING EQUATIONS
Poisson's equation and current continuity equations are governing equations for bipolar ionized field [25]. ∂ρ where ϕ(t) is the total electric potential distribution (V). ρ + (t) and ρ − (t) represent the absolute value of positive and negative space charge density distribution respectively (C/m 3 ). 0 is the vacuum permittivity, whose value is 8.854×10 −12 F/m. J + (t) and J − (t) are current density vectors for positive ions and negative ions respectively (A/m 2 ). R represents the coefficient of recombination and it is approximated to 2.2 × 10 −12 m 3 /s [25]. e is the charge of one electron and its value is 1.602 × 10 −19 C. The current density vectors used above are defined by (4) and (5).
where v + (t) and v − (t) represent the velocity of the positive ions and velocity of the negative ions, respectively (m/s). VOLUME 8, 2020 E(t) represents the electric field distribution (V/m). Note that without other external forces or effects, the negative charges move in the opposite direction of the electric field. W (t) represents the wind velocity distribution of the discussed domain (m/s) and it is a vector. k + and k − are the positive ion mobility 1.4 × 10 −4 m 2 /(V · s) and negative ion mobility 1.8 × 10−4m 2 /(V · s) [25]. The origin and the derivation of the two current continuity equations (2) and (3) are explained in [19]. One helpful note is that the left-hand side of the current continuity equation for negative charges describes the change rate of negative charges while the divergence of J − (t) represents the change of positive charge.
At last, the total current density vector J(t) will be the sum of two current density vectors as shown below.
C. BOUNDARY CONDITIONS Boundary conditions are necessary prerequisites to solve the three coupled governing equations with parameters ϕ(t), ρ + (t) and ρ − (t). We will discuss the definitions of two commonly-used boundary conditions in this section while how to apply them into the proposed algorithm will be introduced in Section III-C. 1) Dirichlet boundary condition is the most explicit boundary condition, which provides the values at the boundary. In this paper, Dirichlet boundary condition may occur on the ground boundary, conductor surfaces of AC and DC lines and also other domain boundaries. Assume that the electrical potential at these boundaries is V 0 and subscript Dirichlet boundary condition with D, the Dirichlet boundary condition can be expressed using the following equation: 2) Neumann boundary condition describes the normal gradient of the unknown functions at the boundary. For this ionized field problem, opposite of the gradient of ϕ(t) is E(t) and so the Neumann boundary condition will provide the value of the normal component of E(t).
Assume the normal component of E(t) is E 0 (E 0 is a known constant) and subscript Neumann boundary conditions with , we will have the equation below.

III. MASSIVELY PARALLEL SIMULATION VIA NODAL DOMAIN DECOMPOSITION
The NDD and upwind NCC methodologies make it possible to fully utilize the computation capability of GPUs. In this section, the details of deploying NDD and applying upwind NCC on GPUs are introduced, followed by the explanation of applying the boundary conditions for solving the Poisson's and current continuity equations. The main idea of this algorithm is to decompose the problem into as tiny sub-domains as possible, so as to maximize the computing efficiency using GPUs' parallel computation ability. To begin, the mesh for the problem is generated by COMSOL Multiphysics TM . Then, the charge density distribution in the domain, including the charge density on the conductor surfaces and the charge density in the air, is initialized. Because of the coupling of governing equations, (1), (2), and (3) are solved alternatingly and iteratively. For each time step, the boundary conditions are applied before using NDD to solve (1) and therefore obtain ϕ(t). When solving Poisson's equation, massive amount of cores in GPUs are executing simultaneously to solve independent sub-domains. By using GPUs, the computation time for a 8184-node FEM problem can be decreased to 30 ms. The upwind difference method is applied to calculate the gradient of positive charge density and gradient of negative charge density. With the known ϕ(t) and gradients of positive and negative charge density, (2) and (3) can be solved simultaneously and the positive and negative charge density distribution for this time step can be calculated. To satisfy Kaptzov's condition, the charge density on conductor surfaces are updated according to E onset and the maximum magnitude of E in the domain. The flow chart of the proposed algorithm is shown in Fig. 1.

A. NODAL DOMAIN DECOMPOSITION
Domain decomposition is an efficient way to solve large-scale problems compared to solving these problems using multiple computational cores in a global system [20], [24]. Following this idea, the nodal domain decomposition (NDD) was proposed in [13]. In NDD, there is only one unknown in each subdomain and each sub-domain can be solved independently and therefore NDD can be easily deployed in GPUs.
The Poisson's equation (1) in Section II-B can be solved by Galerkin FEM, which follows the steps: domain decomposition using FE, elemental formulation by appropriate shape function, global system matrix formation by assembling, the imposition of boundary conditions and solution of the linear system of equations. Same as Galerkin FEM, domain decomposition and elemental formulation are required by NDD. However, NDD avoids the assembly in Galerkin FEM and therefore no huge matrix is involved in NDD and it is a matrix-free method. The system of equations for each element in NDD is generated by integrating the product of the residual and weight function over each element and setting the integral to zero. The elemental equations when linear interpolation (assume ϕ distribution in element (e) can be expressed by ϕ (e) = N 1 ϕ 1,(e) + N 2 ϕ 2,(e) + N 3 ϕ 3,(e) where N 1 , N 2 , and N 3 are interpolation functions) is used are shown in (9).
where K (e) is the coefficient matrix and the calculation of matrix entries for the linear element (triangle) is shown in (10). (e) represents the area of the element (e).
Note that equation (9) can also be expressed by a system of equations although matrix form is used here to keep neat. NDD is a matrix-free method because matrix manipulation or calculation is not involved during the solution process. According to [13], one linear equation including the to-beupdated inner node and all its neighboring nodes can be obtained (in Fig. 2). The values for all the neighboring nodes can be considered as known boundary conditions when updating the inner node. The linear equation for each node is completely independent, which can be solved by GPUs at the same time. Massive parallelism can be achieved by NDD to solve (1). It takes 30 ms to solve Poisson's equation (1) for a FEM problem with 8184 nodes by NDD.

B. UPWIND NODAL CHARGE CONSERVATION
The nodal charge conservation (NCC) guarantees charge conservation law at each finite element node. Besides, each node can be projected to one CUDA core in GPUs to enable solving current continuity equations at each node simultaneously. As a result, the computation efficiency can be improved.
Since the central difference-based iterative scheme usually brings undesired instability [6], the upwind method is utilized. The scenario of the ions migrating in an electric field is quite similar to the situation where a leaf is blown by the wind. The upwind can be defined as the direction the leaf is coming from and so the wind forces the leaf to move downwind. Similarly, in the electric field, we can define the upwind for positive charges or negative charges as the direction the positive or negative charges are coming from.  In other words, the upwind is the opposite direction of the charges' velocities.
Derived from the upwind concept, the gradient of positive or negative charge density at one node is mainly determined by the gradient of the positive or negative charge density in the upwind direction. Then, for one specific node, we can define the element in the upwind direction for positive charges (negative charges) as the upwind element for positive charges (negative charges) (Fig. 3).
With ∇ · W (t) = 0 and equations (1) -(5), nodal charge conservation equations (11) and (12) are obtained by the following three steps. One subscript i can be added to each parameter in (2) and (3) to represent the corresponding parameter for each node and the validity of the equations remains. Then, the first-order discretized forms of the lefthand side of (2) and (3) are used to approximate the partial VOLUME 8, 2020 differentials. At last, the upwind difference method is applied to approximate the gradients of the positive and negative charge densities. Subscript ue represents the upwind element for the corresponding node's positive charge or negative charge [12].
In (11) and (12), the first term of the right-hand sides consists of the gradient of the positive charges or gradient of negative charges and this is where the upwind concept is applied. Note that the upwind element for positive charges or negative charges are determined by the direction of the charges' velocities, which means knowing the charges' velocities is the prerequisite to figuring out the upwind element. In (4) and (5), it is obvious that the charges' velocity has two components: the velocity caused by the electric field and the velocity provided by the wind. The wind vector will be a predefined wind distribution. However, the electric field at each node may be discontinuous because of the fact that electric field E(t) is the gradient of ϕ(t) which is numerically estimated by the linear interpolation.
To provide a reasonable electric field E(t) value at each node, the different averaging scheme has been examined, including the angle-weighted scheme, area-weighted scheme, and unweighted-average scheme. It turns out that the unweighted-average scheme gives the best match with the COMSOL Multiphysics TM result.
From (11) and (12), the computation of charge density for all the nodes in the domain is independent, which makes upwind NCC achievable. With upwind NCC, the oscillation is removed and the parallelism is guaranteed. Solving current continuity equations for each node at the same time significantly improves the computation efficiency.

C. APPLYING BOUNDARY CONDITIONS
This algorithm is suitable to solve the ionized field with both Dirichlet conditions and Neumann conditions. The  Dirichlet conditions can be applied by assigning the known desired value to ϕ at the corresponding nodes but not updating the values at these nodes in the computation process.
However, Neumann conditions cannot be applied directly. Generally, Neumann conditions are described as in (8) with the projection of ϕ in the normal direction. To apply these conditions to the proposed algorithm, the partial differential equation (8) is expanded as a product of two vectors in (13): one of the vectors is the gradient of parameter ϕ at the boundary and the other is the outward unit normal vector of the boundary. On the other hand, according to the Galerkin approach to method of weighted residual, (14) can be obtained from (1) for each element [18]. In (14), it is obvious that the second term on the right-hand side is related to the boundary conditions. When only Dirichlet conditions are applied on the boundary, the integral of the second term is always zero. However, this integral is nonzero when the Neumann condition is applied and the integrated element (e) is a boundary element. The boundary element here means the element who has at least one side coinciding with boundary (the global boundary for the entire problem domain). Assuming boundary condition in (13) is applied and a boundary element (e 1 ) has only one side between local node 1 and local node 2 coinciding with global boundary , then the integral relevant to boundary conditions in (14) can be simplified as (15).
whereâ x andâ y are a set of base vectors in the problem domain and they are in the direction of x-axis and y-axis respectively.â n is the outward vector with unit length which is normal to the boundary andâ n = n xâx + n yây .
where N i and N j are the interpolation functions. Since linear interpolation is used, there are three interpolation functions for each element. (e) represents the boundary of element (e). n x and n y are the projections of vectorâ n normal to the boundary in x and y directions.
where l 12 is the length of the side between local node 1 and local node 2 for element (e 1 ). (e 1 ) is the boundary of element (e 1 ).

IV. CASE STUDY AND RESULTS
The hybrid ionized field around HVDC and HVAC transmission lines whose transmission tower is configured as in Fig. 4 is analyzed. The three-phase 380 kV AC lines are on the left of the tower ordered by phase A, phase B, and phase C from top level downwards. ABC phase sequence is applied and therefore time-varying AC voltages are v A = 310.269cos(ωt) kV , v B = 310.269cos(ωt − 120 • ) kV and v C = 310.269cos(ωt + 120 • ) kV . The bipolar 500 kV DC lines are on the right of the tower ordered by positive polarity, VOLUME 8, 2020 neutral polarity and negative polarity from top level downwards. The voltages of these polarities are V P = +500 kV , V Neutral = 0 kV and V N = −500 kV . Because of the gravity, conductors will sag between transmission towers, which means the distance between the conductors and ground may change -specifically shorten -in the pathway of transmission lines. Short distance signifies a stronger electric field when the stressed voltage is the same. In this perspective, transmission lines where the largest sag happens have the most severe influence on the environment and thus analyzing the ionized field in this location will be more valuable and instructive. According to the Occupational Safety and Health Administration (OSHA), the minimum clearance distance for transmission lines up to 500 kV is 35 feet (around 10.668 m) [16]. In this case, assuming that 10.8 m above the ground is the height of the lowest-level transmission lines with the largest sag. Then, the heights of the three levels of transmission lines after sag are 10.8 m, 18.8 m and 26.8 m (as shown in Fig. 5). The area of 100 m * 80 m (W * H) is chosen as problem domain in Fig. 5 [4].
To solve the hybrid ionized field in the problem domain, the mesh grid was generated by COMSOL Multiphysics TM v5.4. By importing the mesh information into the proposed algorithm which was deployed into GPUs with the CUDA Toolkit 9.1 [1], the case can be solved. A parallel workstation composed of multicore CPUs and many-core GPUs were utilized to solve this case. The CPUs are dual Intel Xeon E5-2698 v4 CPUs@2.2 GHz, 20 cores each with 128 GB RAM. The GPU is NVIDIA Tesla V100-PCIE-16 GB with 5120 CUDA cores [2]. In this case, the time-step t is set to 5 µs in order to assure the stability of the algorithm. The stability condition of this algorithm is that the time-step must not exceed the traveling time for both positive and negative ions in one element [12]. In addition, (16) is applied to update the charge density on the conductor surfaces [22] and E + onset (onset electric field on the positive conductor) and E − onset (onset electric field on the negative conductor) are calculated based on (17) [26]. The parameter R in the equations is the radius of the conductor.
To measure the accuracy and efficiency of the proposed algorithm, the performance of solving Poisson's equation with known charge density distribution by COMSOL Multiphysics TM MUMPS solver and the proposed algorithm is compared. The result comparison of ϕ distribution at each node is demonstrated in Fig. 6 and the mean relative error for all nodes is 0.07%. It turns out that the proposed algorithm is able to provide sufficient accuracy when solving Poisson's equation. In view of speed, the proposed algorithm is 17 times faster than COMSOL Multiphysics TM MUMPS solver when using 40 cores. It is worth mentioning that COMSOL Multiphysics TM MUMPS solver is state-of-the-art [17]. For this 8184-node FE problem, the proposed algorithm achieves to take 30 ms to solve Poisson's equation because of the massive parallelism. Fig. 7 shows the magnitude of the electric field at one node on the DC positive polarity conductor surface and that of another node on DC negative polarity conductor surface with respect to time. Obviously, Kaptzov's condition is guaranteed because either of the magnitudes of E does not exceed the corresponding onset value. Fig. 7 also illustrates the multiple occurrences of corona on the DC conductors. Each time the magnitude of E reaches the onset value, the corona will happen and charges will be generated around the conductors. These charges have a strong effect on weakening the E around conductors because of the proximity and repellency of like charges and therefore a steep drop of E occurs. Due to the diffusion effect, the charges around the conductors go farther away from the conductors and thus the effect of weakening E becomes weaker, then the E on conductors increases slowly. The corona will not happen until the E reaches the onset value for the next time. Fig. 8 demonstrates the same process. We can see two yellow rings in almost all the subplots. One is close to the conductors and the other is far away from conductors. The second yellow ring is due to the second occurrence of the corona. Besides, Fig. 7 demonstrates that the E strength at DC conductor surfaces is composed of DC component and AC component, which is consistent with the statement in [27].
To test the condition with the existence of wind, vector W in equations (4) and (5) is set to 5 m/s and 10 m/s upward to get the charge density distribution in Fig. 8. Refer to (4) and (5), the wind vector will have the same effect on the speed of the charges no matter it is positive or negative. When the wind with 5 m/s upward speed is applied, the charge distributions for both positive and negative charges are not symmetric anymore. The speed of charges moving upward is enhanced while the speed of charges moving downward is restrained. This effect is more significant when the speed increases to 10 m/s. Based on (2) -(5), if the speed of charges changes after applying wind, the current density vectors (J + and J − ) will vary and further the space charge density distributions will be influenced as well.

V. CONCLUSION
In this paper, the massively parallel processing algorithm consist of fine-grained NDD and upwind NCC successfully improves the efficiency of solving the time-varying hybrid ionized field which is computationally burdened with coupling of the Poisson's equation and current continuity equations and repetitive FEM computation. Time-varying distributions of ϕ(t), positive charge density ρ + (t) and negative charge density ρ − (t) can be obtained at each time-step to describe the hybrid ionized field around transmission lines under different wind conditions. With the intrinsic attribute of NDD and NCC to calculate ϕ and ρ for each node independently, the computing power of GPUs can be fully explored to accelerate computation speed. It turns out that the Poisson's equation with 8184-node in the case study can be solved in 30 ms, which is 17 times faster than solver in COMSOL Multiphysics TM , implying the execution time for 10 5 timesteps is reduced to 50 mins instead of 14 hours. Besides, the accuracy of the proposed method is also validated against the commercial software with a mean relative error of 0.07%. The wind will impel or suppress the movement of charges when the wind is blowing to the same direction of the movement of charges or wind is blowing the opposite to the movement of charges.
In the future, this algorithm can be tuned to model other configurations in power systems. In addition, this algorithm can be applied to improve the efficiency of solving similar equations in semiconductor simulation field [9]. Furthermore, the transient finite-element analysis of semiconductors can be solved by proposed algorithm efficiently.