The following premises are adopted for simplicity.

A) All parameters are time-invariant 2-D (only have tow components on the cross-sections of the lines); only uniform corona discharge is considered mainly because real transmission lines are arranged in parallel. Thus, it is assumed that a two-dimensional problem of the corona field is solved on the premise of neglecting poles and towers, conductor sag, and wire end effect.

For the bipolar direct-current ionized field, the equations to be solved include:

Poisson Equation: ∇ 2 Φ = − ρ + − ρ − ε 0 . (1)

Ion Current Equation: j + = ρ + − k + ∇ Φ + W , (2) j − = ρ − − k − ∇ Φ − W , (3)

Current Continuity Equation: ∇ ⋅ j + = − R ρ + ρ − / e , (4) ∇ ⋅ j − = R ρ + ρ − / e , (5)

Φ is the potential; ε 0 is the permittivity of vacuum; j si the current density; k is the migration rate; W is the wind velocity; ρ is the charge density; R refers to the cationic and anionic recombination rates; e is the electron charge.

Calculation processes are shown in Figure 2. First, assumming the charge density on boundary points is constant; then, the electric field and density of space charge are iteratively solved until the steady-state solution under the boundary condition of this charge density is obtained; finally, based on Kaptzov hypothesis, the assumed boundary condition of the surface charge density on the wire are modified, and the space charge density spreading is calculated repeatedly until Kaptzov hypothesis is met.

The electric potential’s boundary condition is expressed as follows.

On the ground and on the Earth wire: Φ = 0 . (6)

On the transmission lines: Φ = V i , i = 1 , ⋯ , s . (7)

To determine the ionized field, the boundary condition of the electric field or electric charge shall also be in place. According to Takuma, supposing that the constant surface charge density on the wire is used as the boundary conditions of the charge, then we get the following equations.

Near the positive wire: ρ + = ρ i , i = 1 , ⋯ , s + . (8)

Near the negative wire: ρ − = ρ j , j = s + + 1 , ⋯ , s + + s − . (9)

In this way, electric charges on the surface of the wire shall be deduced through experiments or empirical formulas before calculations are made, thus reducing the accessibility or reliability of the calculations.

In this study, the field strength remains unchanged at the corona field strength on the surface of the wire.

Near the positive wire: E + = E c i + , i = 1 , ⋯ , s + . (10)

Near the negative wire: E − = E c j − , j = s + + 1 , ⋯ , s + + s − . (11) where E c + and E c − are the corona onset field on the conductors, respectively.

The initial corona field strength of AC lines is calculated by Peek via a great number of experiments. If it is believed that corona onset field of both DC and AC power lines have the same peaks, so Peek’s formula can be transformed into the following DC form. E c = 30 m 1 + 0.301 r k V / c m , (12) where m refers to the surface roughness coefficient, which is generally 0.47–0.5 for the stranded wire; r represents the radius of sub-conducting wire, which is described in the unit of cm.

Eq. 12 presents the results under standard atmospheric pressure. In the case of high-altitude areas with lower atmospheric pressure, the correction coefficient of the atmosphere shall be considered, namely, δ = 0.386 P / 273 + t , where P is the air pressure (mmHg); t is the temperature (°C). E c = 30 m δ 2 3 1 + 0.301 r k V / c m . (13)

There are three independent unknown variables in the bipolar direct-current ionized field. Here, Φ , ρ + and ρ − are selected for simplifying and combining Eqs. 14, 15: V + g r a d ρ + = − k + ε 0 ρ + 2 + k + ε 0 − R e ρ + ρ − , (14) V − g r a d ρ − = − k − ε 0 ρ − 2 + k − ε 0 − R e ρ + ρ − , (15) V + = − k + g r a d Φ + W , V − = k − g r a d Φ + W ,

To calculate the ionized field, by adopting the boundary condition, Poisson Equation and the above two equations are adopted for computing the unknown variables, namely, potential, and space charge density. Iterative solving can be used for the calculation. First, the charge density distribution is fixed, a space-based electric field under certain charge density distribution is obtained and the charge density is calculated. This process is repeated until the field and charge density distribution reach a stable state, as shown in Figure 3.

When solving the electric potential Φ , due to the large gradient of surface potential on the wire, there is a huge error if the FEM is directly used. Here, the electric potential Φ of the total field is divided into the nominal and charge field to solve them respectively; Φ = ϕ + φ . Therefore, Poisson’s equation is transformed into: ∇ 2 ϕ = 0 , (16) ∇ 2 φ = ρ − − ρ + / ε 0 , (17) where ϕ has the same boundary condition with Φ , but φ satisfies φ = 0 on the ground and at the artificial boundary.

In terms of nominal field, multiple images or charge simulation methods are generally adopted for rapidly and accurately solving Poisson Equation. In this study, finite element meshes are divided into triangles, and the same unit partition is also applied to the subsequent charge solving process.

To solve the electric charge, the gradient of node potential shall be obtained, namely E = − g r a d ϕ − g r a d φ , wherein the electrical potential of any point in the nominal field can be directly acquired through the charge simulation method, while the FEM method can only be used for getting the average differential electric field within the triangle units in the charge field. Therefore, the electric field information of the node can be obtained by averaging the electric field information of the node’s surrounding units in accordance with a certain weight. The weight is distributed based on the angle used or the unit area. The calculation results suggest that the angle-based distribution method (as indicated in Figure 4) offers better convergence. In Figure 4, i, j, k, m is the number of computing grid nodes.

To avoid numerical oscillation, upwind discrete equations are adopted while differentially solving the charge density equations. Physically, it can be considered that charge density flows under the electric field effect, and the value at a certain point i is only under the influence of upstream units, as detailed below.

Taking Figure 5 as an example, V i + and V i − are the mobilities of cationic and anionic charges at the node i , respectively. After first-order upwind difference is conducted on positive charges, it is considered that the potential gradient at the point i refers to that within the triangular unit Δ i j k . The first-order difference equation of the triangle unit can be obtained: u ¯ x , y = α 1 + α 2 x + α 3 y . (18)

The coordinates and potentials of the three points are written into: U = u 1 u 2 u 3 , (19) x e = 1 x 1 y 1 1 x 2 y 2 1 x 3 y 3 . (20)

The difference shape function within the unit can be expressed as: α = x e − 1 U . (21)

Thus, the potential gradient is obtained: g r a d U = α 2 , α 3 , (22) − k + ε 0 ρ i + + k + ε 0 − R e ρ i − ρ i + = V i + g r a d ρ i + , (23)

Supposing that x e − 1 = a 1 a 2 a 3 b 1 b 2 b 3 c 1 c 2 c 3 . (24)

The calculation equation of positive charges within the unit Δ i j k can be discretized as follows: V i + g r a d ρ i + = V i + b 1 , c 1 ρ i + + V i + b 2 , c 2 ρ j + + V i + b 3 , c 3 ρ k + . (25)

Finally, the discretized positive charge calculation formula is obtained: − k + ε 0 ρ i + + k + ε 0 − R e ρ i − − V i + b 1 , c 1 ρ i + = V i + b 2 , c 2 ρ j + + V i + b 3 , c 3 ρ k + . (26)

Likewise, the following expression can be deduced for positive charge density: − k − ε 0 ρ i − + k − ε 0 − R e ρ i + − V i − b 1 ′ , c 1 ′ ρ i − = V i − b 2 ′ , c 2 ′ ρ l − + V i + b 3 ′ , c 3 ′ ρ m − . (27)

According to the node information within the unit Δ i j k in the above expression, the updated charge-density value for the node i can be calculated. It is supposed that the surface charge density becomes known in the program, and the charge distribution within the entire space can be solved outwards one by one. Takuma suggested that the binary second-order equations with positive and negative charge densities as the variables are directly solved, and the following constraints are used for determining the reasonableness of the equation root:

A) If ρ j + ≥ 0 , ρ k + ≥ 0 , and ρ i − ≥ 0 , 0 ≤ ρ i + ≤ max ρ j + , ρ k + , ρ i − ;

From the physical perspective, it can be explained that the charge density within a space should not be greater than that on the boundary of the wire. However, direct solving of equations heavily depends on the selection of the initial value of space charge, whose convergence cannot be guaranteed. In this study, the equations are differentiated into the first-order form, and the following equations are iteratively solved: ρ i , n + 1 + = V i + b 2 , c 2 ρ j , n + 1 + + V i + b 3 , c 3 ρ k , n + 1 + − k + ε 0 ρ i , n + + k + ε 0 − R e ρ i , n − − V i + b 1 , c 1 , (28) ρ i , n + 1 − = V i − b 2 ′ , c 2 ′ ρ l , n + 1 − + V i + b 3 ′ , c 3 ′ ρ m , n + 1 − − k − ε 0 ρ i , n − + k − ε 0 − R e ρ i , n + − V i − b 1 ′ , c 1 ′ . (29)

Estimation of Initial Value of Charge Density on the Boundary of the Wire.

In this study, the Kaptzov hypothesis is adopted as the boundary condition. Before the calculation is made, assuming that there are some charges distributed, and the electric field within the space is calculated. Then, the surface charge density distribution on the wire is updated based on the difference between the field strength and corona onset field.

The initial value of surface charge density on the wire is assumed as follows; the expression can be obtained near the positive wire: ρ + i , 1 = ρ e ⁡ cos π − θ i / 2 , i = 1,2,3 , ⋯ , M , (30) ρ − i , 1 = 0.1 ρ + i , 1 , i = 1,2,3 , ⋯ , M , (31) where ρ e = E y D E c r i t R 8 ε 0 Φ c Φ − Φ c / D 2 Φ 5 − 4 Φ c / Φ , (32)

It is the same with the expression near the negative wire. E y refers to the field strength of the point at the height H on the axis y. Φ c refers to the corona voltage: V 0 = m r 2 E c ln 2 H r . (33)

Takuma takes the initial value of charge density of each mesh on the surface of the wire as ρ e , but for real bundle conductor lines, the electric field strength is very unevenly distributed on the wire surface, and corona discharge intensity also differs. If the same initial value of charge density is selected, it is not conducive to computing stability and efficiency. Therefore, the discrepancy between surface field and corona onset strength is chosen as the benchmark to estimate the initial value of space charge.

First, the maximum nominal surface field strength E max is calculated. The initial value of space charge for point i on the wire surface is estimated.

For positive wire: E max ≥ E c : ρ + i , 1 = ρ e E i − E c / E max − E c , (34) ρ − i , 1 = 0.1 ρ + i , 1 , (35) E max < E c : ρ + i , 1 = ρ − i , 1 = 0 . (36)

For negative wire: E max ≥ E c : ρ − i , 1 = ρ e E i − E c / E max − E c , (37) ρ + i , 1 = 0.1 ρ − i , 1 , (38) E max < E c : ρ + i , 1 = ρ − i , 1 = 0 , (39)

For charge correction on the wire surface, the difference between surface field and onset field of the wire is adopted, which resembles the secant method, as described below: ρ + i , n + 1 = ρ + i , n 1 + μ E i − E c / E i + E c , i = 1,2,3 , ⋯ , M , (40) where E i and E c refer to the surface field strength and onset field strength; μ is the updated accelerating factor, which is equal to 2 in this study.

The initial value of space charge density except for those on the boundary is set as 0.

Taking the ± 800 kV Hami-Zhengzhou power line at the high altitude (with an altitude of 2,000 m) as an example, the wires are horizontally arranged in a bipolar way, and insulator strings adopt a V-shaped layout; the straight-line tower has a model of Z27102A1, and 6*JL/G3A-1000/45 steel-cored aluminum strand is used as the wire; JLHA1/LB1-180/30 aluminum alloy conductor aluminum steel reinforced strand is adopted as the ground wire.

Major parameters of wire and ground wire are shown in Tables 1, 2.

Other parameter settings for the program are shown in Table 3.

When typical ± 800 kV power lines in the high-altitude region have different heights above ground, the ground field distribution are shown in Figures 6, 7.

According to relevant national standards, for ± 800 kV UHVDC overhead transmission lines, in case of maximum sag, the minimum distance between the wire and the ground is 21 m when the lines go through the residential areas, and non-distorted total electric field on the ground shall not exceed 15 kV/m. It can be seen from Figure 7 that the total electric field on the ground increases with the wire’s decreasing height above ground. 6*JL/G3A-1000/45 steel-cored aluminum strand is adopted; when the wire height reaches 21 m, the electric field on the ground is 10.2 kV/m, which is much lower than the limit set in the national standards.