Fast Corona Discharge Assessment Using FDM integrated With Full Multigrid Method in HVDC Transmission Lines Considering Wind Impact

A novel approach for solving the monopolar corona in high voltage direct current (HVDC) transmission line systems is proposed by the ﬁnite difference method (FDM) and a full multigrid method (FMG). Speciﬁcally, the FMG is implemented as a fast solver with respect to existing iterative solutions for the FDM to solve the Poisson equation, particularly on ﬁne grids. The advantage features of the proposed approach are that it avoids the hypothesis of a constant electric ﬁeld around the conductor’s surface. Further, it considers the inﬂuence of space charge on both the magnitude and the direction of the electrical ﬁeld. The proposed approach is employed for computing the electric ﬁeld and current density on the ground plane with and without wind effect. Considering the impact of wind in the present study, the ﬁndings conﬁrm that both corona current density and electric ﬁeld on the ground plane are inﬂuenced by the transverse wave. Eventually, the effect of changing the wind speed on the electric ﬁeld proﬁles and the current density is deeply studied in HVDC transmission line systems. To prove the efﬁcacy of the proposed approach, it is compared with previous experimental results where a better agreement is reached rather than other numerical techniques.


I. INTRODUCTION
With the huge increase in power transmitted over long distances, high voltage direct current (HVDC) transmission has become emulative, and many HVDC transmission lines have been built around the world [1], [2].HVDC transmission uses the ground as a return path, requiring simpler towers and therefore lower net costs.However, HVDC transmission systems operate above the corona onset voltage and may cause audible noise, radio interference, power loss and may lead to deterioration of the insulation systems [3].In addition, corona might cause biological and environmental effects.This phenomenon has formed a growing concern regarding The associate editor coordinating the review of this manuscript and approving it for publication was Sonia F. Pinto .
the study of corona discharges on HVDC transmission lines.
In HVDC transmission systems, corona generates space charges that fill the region between the electrodes.However, in high voltage alternating current (HVAC) transmission systems, space charges produced by corona are limited to the proximity of the conductor as the electrical field is periodically reversed.Corona loss depends mainly on the conductor's surface field in the case of HVAC transmission systems, while in the case of HVDC transmission systems, it depends more on the geometric parameters of the transmission line [4].Due to these reasons, the study of the electric field and the current density associated with corona discharge, in HVDC transmission systems, is of great interest to researchers.Generally, the corona effect should be considered in the design of HVDC transmission systems.
The corona problem in the vicinity of the HVDC transmission systems is of a non-linear nature, and so the exact solution of the corona-associated ionized field is extremely challenging.In this regard, several attempts were made to develop numerical techniques.The main purpose of these techniques is to solve Poisson and current continuity equations while obtaining the electric field and the current density associated with corona on monopolar HVDC transmission systems [5].Common numerical techniques for solving the corona problem are: 1) the finite element method (FEM) [6]- [14], 2) the charge simulation method (CSM) [15]- [18], and 3) the boundary element method (BEM) [19], [20].
Sarma and Janischewskyj [21], [22] developed a computer model for the analysis of HVDC transmission lines.In particular, they proposed an iterative solution for the electric field and the current density on the ground surface.The obtained voltage-current characteristics were reasonably consistent with the experimental data where their analysis implemented Deutsch assumption which mentions that the space charge influence only the magnitude of the electric field, but not its direction The validity of the Deutsch assumption has remained an existing issue to be studied by the researchers.Takuma et al. [9] used the FEM to solve the ionized field problem for the mono-polar corona, assuming that the space charge density is constant around the conductor surface.They abandon the assumption that the electric field remains constant at its onset value on the conductor surface.Nevertheless, the choice of the constant charge density around the conductor surface as a boundary condition is not practically realistic, which puts some empiricism on their analysis by the FEM.Abdel-Salam et al. used the FEM in solving the monopolar ionized field eliminating both Kaptzov and Deutsch assumptions [7].Davis and Hoburg suggested a numerical solution based on the FEM combined with the method of characteristics using Kaptzov's assumption in their analysis [6].Abdel-Salam and Al-Hamouz came up with an iterative FEM for the solution of the ionized field problem in monopolar transmission lines taking into consideration the ion diffusion [10].
The finite difference method (FDM) is a numerical technique that provides approximate equations by neglecting higher-order terms, which may cause an excessive truncation error.This problem could be addressed by means of finer grids.However, this will increase the computational time relative to the other conventional iterative methods.The convergence rate of the most common iterative approaches, such as the Jacobian approach, the Gauss-Seidel approach, and the successive over-relaxation (SOR) approach, is low as the grid becomes finer.In such cases, the need for another efficient iterative method is necessary to implement finer grids without suffering from high computational time.This problem is addressed by the full multigrid method (FMG) and the convergence rate of the traditional iterative approaches is increased [23]- [25].Generally, the concept of FMG is based on converting the fine mesh into a coarser mesh.On the coarser mesh, one of the conventional iterative techniques should be implemented to handle the coarser mesh.After that, the solution obtained is utilized by grid transfer operators as an enhanced initial estimation for the fine mesh [23].Recently, the FMG is applied as a new iterative solver for the FDM in resolving the corona phenomenon in the wire duct precipitators, and it improves the computational time of the FDM compared to the classical iterative methods, especially when using fine meshes [26].The FDM with the FMG showed promising results for solving the corona problem in the wire duct precipitators.However, its usage for solving the monopolar corona in HVDC transmission line systems is not yet investigated, which is covered in this work.Further, most of the previous studies are missing to analyze the effect of changing the wind on the electric field profiles and the current density in HVDC transmission lines by the FDM with the FMG solver.
To cover the gap in the literature, this study is aiming to propose a novel approach for solving Poisson and the current continuity equation of the monopolar corona in HVDC transmission line considering wind effects.This approach is accomplished by combining the FDM with the FMG to operate on the finer mesh without suffering from high computing time.The proposed analysis avoids both Kaptzov's and Deutsch assumptions.In addition, the present work analyzes the nature of the corona problem in the presence of wind and studies its effect on the electric field and the current density values above the ground plane.Finally, the effect of changing the wind speed on the ionized field and the current density distribution is presented.The proposed numerical technique is compared with previous published numerical techniques and experimental work including no wind and wind conditions.
The major contributions of the paper can be listed as follow: • The other published numerical techniques have not dealt with finer computational domains, as this will increase the computational time.
• Unlike the other published techniques in the field of HVDC systems, the present technique proposes the FMG for the first time to overcome excessive computational time especially on solving the Poisson equation in the ionized field problem on dealing with fine computational domains.
• For the first time in solving the Poisson equation, the FMG is introduced as an iterative solver for the finite difference equations, and it gave the FDM the advantage of working on fine computational domains without suffering from the high computational time.
• Working on finer grids gives a more accurate picture for field and current density distribution which agreed better with the experimental results than the other published techniques and this has been proved in the comparisons.
• Using the finer grids for a certain design of HVDC systems would help a lot in getting an accurate picture of the electric field and density distributions and can be  more confident about the predictions of the numerical outcomes, without the need of the costly experiments.So, it can be an effective tool for engineers in the early design stage of HVDC mono-polar systems.

II. PROPOSED CALCULATION METHOD A. CORONA GOVERNING EQUATIONS
The following equations define the discharge of the monopolar DC corona as follows [1]: Poisson's equation; Current continuity equation; where E is the electric field intensity, j is the current density, V is the electric potential, k is the ionic mobility, is the air permittivity, W is the wind velocity, and ρ is the charge density.Combining ( 1)-( 4) results in a non-linear third-order equation [7] leading to significant difficulties in the analytical solution.Therefore, numerical techniques play a significant role in the treatment of the monopolar corona problem.

1) BOUNDARY CONDITIONS
Referring to Fig. 1, the boundary conditions necessary for solving the mono-polar DC corona equations are: 1) The ionized wire surface potential is equivalent to the applied voltage (V ).
2) The potential on the ground is zero.
3) The space charge in the inter-electrode region is unipolar.4) The thickness of the ionization region around the ionized wire is ignored compared to the height of the wire (H ) above the ground.
5) Along the line of symmetry and at the ground plane, the electric field E x equals zero.6) The ion mobility k is constant.7) The charge density at the point A (ρ A ) can be calculated from [26]: in which f represents the factor of the surface roughness, δ represents the factor of air density, R represents the ionized wire radius (cm), j p represents the mean current density on the plate understudy, S x represents the half-length of the grounded plate, and k represents the ion mobility (m 2 /Vs).
From a mathematical point of view, the solution of the monopolar corona problem, defined by equations ( 1) to (4), requires three boundary conditions for any conductor geometry.It is required to provide sufficient estimations of the corona inception field of line conductors [27].In the present study, the assumptions of the electric field constancy at the ionized wire surface at its onset value (Deutsch assumption), and that the space charge impacts only the magnitude, but not the direction of the electric field (Kaptzov assumption) are ignored.Therefore, the present study inserts ( 5) as a boundary condition to give a certain current density on the plate understudy besides that the ionized wire surface potential is equivalent to the applied voltage (V ), and the potential on the ground is zero.

2) SOLUTION OF THE CONTINUITY EQUATION
The FDM is implemented to solve the continuity equation in the case of wind and without wind conditions.Combining ( 2) and (3) should result in: Substituting W by zero in (6) leads to: Considering that the ion mobility is constant and substituting in (7), the current continuity equation will result in: Applying backward finite difference equations, equation (8) will give: So, the charge density at any point ρ (i,j) is computed by solving the following equation, in which E x and E y represent, respectively, the electric field intensity along the X and Y directions.Also, x and y represent, respectively, the incremental distance along with X and Y directions.Taking into consideration, the wind velocity vector in solving the continuity equation by finite difference method, equation ( 6) becomes: Substituting ( 1) in (11) results in: As the wind is divergence less, i.e.
∇.W = 0 (13) So, equation ( 12) ends to, The wind velocity vector is given by, W = w x a x + w y a y (15) In the present study, the proposed method deals only with a wind velocity component in X direction, so, equation ( 14) should give, Using the backward difference formulas, the above equation becomes: So, the charge density at any point ρ (i,j) is found by solving the following equation: Then, the current density on the grounded plate could be computed from Obviously, the presence of wind will shift every point in the grid depending on the magnitude of the wind velocity and electric field intensity in this region.Therefore, the resulted grid will be un-equally spaced, which should satisfy: in which E x and E y represent the electrostatic fields in the horizontal and vertical directions, respectively.

3) SOLUTION OF POISSON EQUATION
The area of the intended domain in the mono-polar configuration is not bounded as shown in Fig. 1.It should therefore be specified by an artificial boundary.This artificial boundary will extend S x in the X -direction and S y in the Y -direction.At these boundaries, the resultant electric potential is assumed to be equal to the electrostatic potential [9].The domain is divided into rectangular elements and the finite difference method is implemented to solve the Poisson equation via the FMG as an iterative technique.
As this paper concerns with studying the effect of wind velocity on the current density and the electric field profiles on the ground plane, the points of the grid will be unequally spaced from each other.So, the finite difference equations should be formulated for unequal spacing in the grid.
Referring to Fig. 2 and using Taylor's expansion in the horizontal direction: by solving (21) and ( 22) the result will be in (23): Similarly, using Taylor's expansion in the vertical direction, Substituting ( 23) and ( 24) in (1), the potential at any point in the grid can be calculated.

VOLUME 8, 2020
For equally spaced points in the grid, (i.e.h e = h w = x and h n = h w = y), and substituting ( 23) and ( 24) in (1) leads to Accordingly, the potential v i,j is determined at each possible node in the grid by applying the FMG solver as an iterative tool.

B. FULL MULTIGRID METHOD
Generally, smoothers, as well as restriction, are the backbone of the multigrid method, besides the prolongation [23].The main aim of smoothers is to reduce high-frequency errors.Among the grid transfer operators, there is a constraint that maps the imperfection in residues from the finer grid to the coarser grid after smoothing the finer grid.Restriction can be implemented through 1) straight injection, 2) half weighting, and 3) full weighting.
The additional grid transfer operator implemented with the FMG is an interpolation process.By employing interpolation, we can calculate the mean value of all present adjacent points according to the amount of each next grid point [23].Obviously, the Poisson equation consumes excessive time in the computational algorithm [26].When the grid is fine, it takes a lot of time to achieve the desired convergence through classical iterative methods for the solution.
The present work, therefore, suggests the multigrid approach as an iterative solver to improve the convergence process.The multigrid approach has different schemes, including 1) the two-grid system, 2) the V-cycle, 3) the W-cycle, and 4) the FMG.In the suggested approach, FMG is applied as an iterative tool for solving the Finite-Difference equations while the Gauss-Seidel technique is executed as a smoother.
The FMG begins with the coarsest grid and then interpolate the result of the finer grid whereas numerous V-cycles are completed.
Considering the Poisson equation ( 1) in the form of: The one V-cycle is implemented as follows [23] 3. Restrict the residual to a coarser grid, r coarse ; where [I h 2h ] is the standard full weighting operator [23].

In the case of the grid is the coarsest one, answer the residual equation (denoted by [B][e] = [r]
) to determine an estimated value for the error [e]; otherwise, return to step 1.

Apply the interpolation technique to correct the coarse
grid to a finer grid : in which [I 2h h ] represents the operator of bilinear interpolation.6. Relax post-smoothing number steps (N2) by the Gauss-Seidel technique smoother on the finest grid to determine alternative initial estimation for the next V cycles.
Several V-cycles are carried out till reaching the finest grid.Next, a deeper V-cycle is completed on the finest level till reaching the preset convergence.Here, it was found that N1 = 3 and N2 = 1, can yield the optimal performance in terms of the problem computational time.
Figure 3 shows the flow chart of the proposed approach that solves the monopolar corona in HVDC transmission line systems by the FDM integrated with the FMG solver.To implement the proposed approach, a Matlab program has been developed to resolve both Poisson and continuity equations.

A. NO IMPOSED WIND RESULTS (W = 0) 1) VARIATION OF ELECTRIC FIELD AND CHARGE DENSITY AROUND THE CONDUCTOR'S PERIPHERY
The scale model of Hara et al. [28], with 0.25 cm wire radius located at a height of 2 m from the ground level, is studied.The ion mobility and the roughness factor are 1.5 × 10 −4 m 2 /Vs and 0.942 respectively as defined by Hara et al. [28].The finite difference method is applied to a fine rectangular grid of (256 × 256) points, with an artificial boundary extended to 2H in the vertical direction (S y ) and 3H in the horizontal direction (S x ).
The proposed method is utilized to evaluate the electric field intensity around the periphery of the ionized wire at a stressed voltage of 300 kV as shown in Fig. 4. The results show that as the applied voltage is more than the corona onset voltage, the electric field around the wire (E w ) drops below the onset field (E on ) in the same manner as was found experimentally in previous work [29], [30].At the point on the ionized wire facing the ground plane, the electric field decreases by 15% below its corresponding onset value.So, the electric field around the ionized wire is not the same at its onset value which affirms that the present study did not assume the constancy of the electric field on the wire surface at its onset value.
Figure 5 demonstrates the variation of the charge density (ρ) around the surface of the ionizing wire when a voltage of 300 kV is applied.The results affirm that the current work does not involve the hypothesis of a constant charge density surrounding the wire periphery.The maximum variation between the maximum and the minimum charge density  around the wire surface is about 0.27%, which make the hypothesis made by Takuma is somewhat reasonable [9].
2) COMPARISON BETWEEN FDM AND EXPERIMENTAL RESULTS OF HARA et al. [28] The present analysis is used to calculate the voltage and current characteristics for Hara's model mentioned in the previous section [28].A strong matching between the proposed method and the experimental data of Hara at various applied   voltages, as shown in Fig. 6, was found, which affirms the accuracy of the proposed approach.
Similarly, the proposed approach is compared with the results of the FEM and the method of characteristics (FEM-MOC) model proposed by Davis and Hamburg [6].The proposed technique fits the experimental data better than [6] technique, especially at higher voltages.Comparison between FDM and previous measured and calculated results at 200 kV [6], [12], [28], [31].(a) Current density distribution along the ground plane, and (b) Distribution of the electric field along the ground plane.
Figures 7 and 8 illustrate the variation of the measured and the calculated values for the current density (J g ) and the electric field (E g ) at the ground plane for Hara's model at 200 kV and 300 kV, respectively.The FDM is compared to previously published work on the model of Hara, and the results fit the experimental data of Hara, especially at the area under the ionized wire on the ground plane, better than the other numerical methods introduced by Davis and Hoburg [6], AL-Hamoz and Zakariya [12], and Adamiak and Kazimierz [31].

3) COMPARISON WITH THE EXPERIMENTAL AND THE CALCULATED RESULTS OF ABDELSALAM et al. [7]
The finite difference approach is also examined against the experimental data of Abdelsalam et al. for a wire of a radius 0.165 cm positioned at a height of 0.4 m over the ground surface with ion mobility of 1.5×10 −4 m 2 /Vs [7].The wire is stressed at voltages between 60 kV and 90 kV.The proposed method computes the voltage and current characteristics and is compared to the FEM results of Abdelsalam et al, and a FIGURE 8. between FDM and previous measured and calculated results at 300 kV [6,28].(a) Current density distribution along the ground plane, and (b) Distribution of the electric field along the ground plane.
better agreement with the experimental results is reached (see Fig. 9).Further, the current density at the ground plane is calculated using the finite difference method for the Abdelsalam experimental model at applied voltages of 80 kV and 90 kV, as given in Figs. 10 and 11, respectively.
Generally, a strong agreement is reached between the finite difference approach and the numerical results of Abdelsalam et al., with a slight exception to the current   that it is higher than the values taken from the V-I characteristics by 14% and 17% for nominal voltages of 80 kV and 90 kV, respectively.e) The ground plane of the experimental setup of Abdelsalam is 0.75 m in length around both sides of the wire, and the current sensing segments were 4.5 m in length equal to the wire length.Consequently, end effects may affect the experimental results.
B. IMPOSED WIND RESULTS (W =0) 1) COMPARISON BETWEEN FDM AND THE EXPERIMENTAL RESULTS OF HARA et al. [28] With considering the wind profile, the FDM is applied to calculate the current density and the electric field on the ground surface for Hara's model, for nominal voltages of 200 kV and 300 kV [28].
Figure 12 shows the variation of the current density and the electric field on the ground surface, for an imposed wind with velocity 8 m/s below the stressed conductor at 200 kV.The wind direction is assumed to blow from west to east.It is VOLUME 8, 2020 found that the highest values of current density and electric field are not at X g = 0 m, as in the no wind condition, but it is shifted to the right in the direction of the wind velocity.The proposed method is compared with the FEM-MOC model introduced by Davis and Hoburg [6], for an ion mobility value of 1.5 × 10 −4 m 2 /Vs.Accordingly, the proposed method fits the experimental data of Hara, better than the method proposed by Davis et al..In addition, the proposed method tests the variation of the current density on the ground surface at ion mobility of 2.1 × 10 −4 m 2 /Vs.It was found that this value of ion mobility gives a better agreement concerning the magnitude and the position of the experimental data of Hara's model.
Also for the same model of Hara et al., the proposed approach is compared with the upstream FEM introduced by Lu et al. [13], and the FDM agrees with the experimental data better than FEM, as shown in Fig. 12.
To prove the validity of the proposed method, the FDM is implemented for the same model at the same wind speed by applying a voltage value of 300 kV.A good agreement is reached with the experimental data and a better fit concerning the magnitude and the position of the peak values for the current density and the electric field at the ground surface compared to the findings reached by Davis et al. [6], as shown in Fig. 13.
It is important to note that the variation of ion mobility affects the distribution of the electric and current density.In the case of Fig. 7, Hara's model is conducted for no wind conditions at ion mobility of 1.5×10 −4 m 2 /Vs, but in the case of Figs. 12 and 13, Hara's model is conducted in the presence of wind.The wind exerts a mechanical force on the ions in HVDC transmission systems, which in turn increases the drift velocity of ions and the ion mobility towards the ground from 1.5 × 10 −4 to 2.1 × 10 −4 m 2 /Vs as compared to experimental results.It was found that ion mobility of 2.1 × 10 −4 m 2 /Vs gives a better agreement with the experimental results by Hara et al., and this value is within range for ion mobilities mentioned in literature.

2) WIND SPEED EFFECT ON THE CURRENT DENSITY DISTRIBUTION
Figure 14 demonstrates the influence of changing the wind speed on the variation of the current density on the surface of the ground for Hara's model at an applied voltage of 300 kV.The change of the wind speed affects significantly the distribution of the current density concerning its magnitude and position of the highest current density resulting in asymmetric distribution.However, the highest current density at the ground surface increases by 80.9%, 86.9%, and 100.5 % compared to no wind conditions, (see Fig. 8.a) for wind speeds of 2, 4, and 8 m/s, respectively.In addition, the location of the maximum current density is shifted to the right as the wind speed increases.It should be observed that the magnitude of the current density in the region located to the left of X g = 0 m, decreases as the wind speed rises.Whereas in the region located to the right of X g = 0 m, the magnitude of the charge density rises with increasing the wind speed.Also, as the wind velocity increases, the dispersal of the current density on the ground plane becomes broader.

3) WIND SPEED EFFECT ON THE ELECTRIC FIELD INTENSITY DISTRIBUTION
For the electric field distribution, the wind speed affects only the of the maximum electric field, but the value of the maximum field intensity is almost constant as shown in Fig. 15 compared to the no wind condition, Fig. 8.b.The value of the electric field at both sides of X g = 0 acts in the same manner as explained in the preceding section for the current density when the wind speed increases.
The above results can be explained that the wind exerts a mechanical force on the ions, which in turn decreases the space charge at the gap central part, enhancing the field above the conductor and gives arise relating to the value of the corona current, while maintaining an approximately constant field intensity on the ground plane, nonetheless of the growth in the current density on the ground plane.
Generally, the FDM with the FMG solver predicts a significant rise in the value of the current density, and a closely persistent electric field on the ground plane when the wind speed increases and these results agreed qualitatively with the theoretical model of Takuma et al. [9], and the previous experimental results [32]- [34].
The power of the FMG is concentrated in the good initial guess for the potential values, besides working on different domain sizes with the aid the transfer operators mentioned in the manuscript.In terms of computational time, the FMG is greatly transcended owing to the timing performance as compared to the classical iterative solvers.For example, using the FMG as an iterative solver, the overall program execution time is 8 and 30 seconds to converge to a prescribed tolerance in the case of the no wind and wind conditions, respectively, on a grid of (256 × 256) points.While in the case of Gauss-Seidel as an iterative technique, it has a poor convergence rate as it takes around 5 hours to converge to a certain tolerance.Therefore, the FMG requires fewer iterations which in turn decreases the computational time.

IV. CONCLUSION
A novel approach for resolving the monopolar corona issue related to the HVDC transmission line configuration using the FDM integrated with the FMG solver has been analyzed.The FMG improves the convergence rate of the FDM, on the finer grids.The overall program execution time is 8 and 30 seconds to converge to a prescribed tolerance in the case of the no wind and wind conditions, respectively, on a grid of (256 × 256) points.The proposed approach is evaluated concerning previous experimental data and a strong matching is achieved, compared to the previous numerical techniques for wind and no wind conditions.The results show that the wind exerts a mechanical force on the ions in HVDC transmission systems.The present study predicts a significant rise in the current density and almost constant electric field on the ground plane as the wind velocity increases and leads to a broader distribution for the current density and the electric field on the ground plane.The future work will be directed to implement the proposed FDM-FMG in bipolar HVDC lines on fine computational domains.

FIGURE 2 .
FIGURE 2. Fractional grid of the calculational domain.

FIGURE 3 .
FIGURE 3. Flow chart of the computational algorithm.

FIGURE 4 .
FIGURE 4. Electric field variation with the onset value around the surface of the wire (θ ).

FIGURE 6 .
FIGURE 6. V-I characteristics of Hara scale model.

FIGURE 7 .
FIGURE 7.Comparison between FDM and previous measured and calculated results at 200 kV[6],[12],[28],[31].(a) Current density distribution along the ground plane, and (b) Distribution of the electric field along the ground plane.

FIGURE 12 .
FIGURE 12.Comparison between FDM and previous measured and calculated results at 200 kV and wind speed of 8 m/s [28].(a) Current density distribution along the ground plane, and (b) Distribution of the electric field along the ground plane.

FIGURE 13 .
FIGURE 13.Comparison between FDM and previous measured and calculated results at 300 kV and wind speed of 8 m/s [28].(a) Current density distribution along the ground plane, and Distribution of the electric field along the ground plane.

FIGURE 14 .
FIGURE 14. Distribution of current density along the ground plane for various wind speeds at 300 kV.

FIGURE 15 .
FIGURE 15.Distribution of electric field along the ground plane for various wind speeds at 300 kV.