Calculation of the electric field distribution in the vicinity of the conductive rod

The article reviews the methods of mathematical modeling of electric fields in the vicinity of conducting rods and presents a method developed for calculating the distribution of the electric field strength and potential in systems with conducting rods. This method allows to use a computational spatial grid with a step proportional not to the radius of the rod, but to its length, which is relevant when the ratio of the rods length to its radius is large. The method is applied to the calculation of rods EF, for which this ratio is of the order of 10–10. The proposed method is based on the finite integration method. At the same time, the nonlinear decrease in the levels of strength and potential when moving away from the rod in directions perpendicular to its axis is taken into account. The difference coefficients at the nodes surrounding the rod were obtained by integrating over the computational grid cell surfaces of expressions describing the strength and potential of the electric field for an elongated conducting ellipsoid under potential. With this representation of the conducting rod, it was possible to achieve the greatest accordance of calculations with the analytical solution. In practice, the application of the presented method allows for a more accurate calculation of the electric field in the vicinity of a conducting rod, which is either under potential, or in a homogeneous electric field, using a computational grid with a step proportional not to the radius of the rod, but to its length. The non-linear character of the decrease in the strength and potential of the electric field near the rod is taken into account using analytical expressions for a conductive ellipsoid under potential. In the area surrounding the rod and above its top, when using a spatial grid step proportionate with the length of the rod, and not with its radius, the relative errors in calculating the strength decreased from 27 % to 3 %. The results of calculating the electric field of a lightning rod are presented in order to analyze the conditions for the occurrence of upward leaders.


Introduction
Solving the problems of modeling rod lightning rods, incomplete insulation breakdown channels, channels of lightning leaders [1][2][3][4][5][6] may be performed by calculation of the electric field (EF) in the vicinity of conducting rods, the length of which exceeds their radius by several orders of magnitude (L/R>10 2 -10 3 ). Another area where the solution of such problems is required is emission devices using arrays of carbon nanofibers [7][8][9].
Determination of the degree of strengthening of the EF tension in the area of the top of the rod is an urgent problem in various areas [10,11].

Analysis of recent research and publications
When calculating the distribution of the electric field in the vicinity of long and thin rods, we encounter a problem associated with the choice of the step of the computational grid. The step should be less than the radius of the rod, however, with a significant ratio between the length and the radius of the rod, the resulting smaller computational grid leads to a significant increase in the order of the equations system describing the system, which makes it impossible to solve it using existing computational tools.
In [12], for three-dimensional modeling of the gas discharge process in order to study the streamer branching processes, an approach was chosen, which consists in utilizing a computational grid divided into several layers. The step in each layer becomes smaller and smaller as it approaches the discharge axis. The solution is based on the finite element method using Newton's iterations to solve a nonlinear system of equations, which ensures the simultaneous convergence of the entire system of equations. Lagrangian-quadratic elements were used for the Poisson equations and boundary conditions. This approach allows you to perform calculations with the required accuracy, but the requirements for RAM and time costs are a significant obstacle to its use.
In the numerical calculation of the EF using finitedifference methods in the immediate vicinity of a thin cylindrical rod having an infinite length [13][14][15][16][17], to determine the coefficients of the difference equations, the law of decay of the EF intensity is used, which is inversely proportional to the distance from its axis. In this case, the step of the computational grid can significantly exceed the radius of the rod. In [18], a rod of finite length is represented by a uniformly charged thread, in [19] by an elongated ellipsoid located in a uniform EF. In both cases, when calculating the area in the vicinity of the rod, the same law of strength variation was applied. When comparing the results of the potential distribution calculation at the distance of the step of the computational grid from the rod by numerical methods and analytical calculation of a uniformly charged thread and an elongated ellipsoid located in a uniform EF [20], a relative error of about 5 % was obtained. The error in calculating the strength exceeds 10 %. When using this approach in solving problems containing a rod under potential, in the absence of an external EF, the error becomes even greater: up to 45 % for the strength calculation of the rod top zone under potential. In order to solve this problem in the region limited by the step of the computational grid from the axis of the rod, when calculating the coefficients of the difference equations at the nodes on the rod and the step from it, it is proposed to use analytical expressions for the EF of an elongated ellipsoid under potential [20].

The aim of thе work
The aim of this work is to develop a refined method for the numerical calculation of the electric field in the vicinity of conducting rods with a large ratio of length to radius. The difference coefficients in the nodes of the computational grid surrounding the rod are determined in a special way. As a result, the used step of the computational grid can be proportional to the length of the rod rather than its radius. This approach increases the accuracy of calculating the EF strength in systems containing conducting rods.

Basic research materials
The calculation used a system ( Fig. 1) with a grounded rod with potential U 0 .

Fig. 1 -Calculation scheme
The computational equations used in the finite volume method [21,22] were obtained by integrating Maxwell's equations over the volumes V of the unit cells that make up the computational domain. The nodes of the computational grid (i, j, k), in which the EF potentials are determined, are located at the interface between the media, and therefore on the axis of the conducting rod ( Fig. 2). With such an arrangement of the computational grid nodes, the conditions at the interfaces between the media are satisfied automatically.

______________________________________________________________
After applying the divergence operation to the equation, the resulting expression must be integrated over the volumes of the unit cells V that make up the computational domain using the Ostrogradsky-Gauss theorem. In the absence of space charges, the second term on the right-hand side of the last equation can be neglected. As a result, an equation will be obtained for determining the electric potential φ: where n -normal to the surface S enclosing the volume V, E n -the projection of the vector strength E normal to the surface S.
For each node (i, j, k) of the computational grid, an equation of the form (1) is drawn up. In this case, each cell of the computational domain is characterized by its specific conductivity γ i j k . Since the length of the rod is several orders of magnitude larger than its radius, the rod in the design model is replaced by a set of nodes located on the axis with indices r (Fig. 2).
The non-linear nature of the decrease in the strength and potential of the electric field in the direction perpendicular to the rod axis is taken into account using the conductivity tensor: where k x , k y , k z -coefficients equal to 1 for all nodes except (i r , j r , k r ), (i r -1, j r , k r ), (i r , j r , k r -1).
The presence of conductivity between the nodes i r , j r , k r along the rod axis is taken into account using the coefficients k y = 10 6 for these nodes.
The coefficients k x and k z of the nodes surrounding the rod and lying on its axis, and the coefficient k y of the node at the top of the rod are determined from analytical expressions for the potential level and the electric field strength depending on the voltage applied to the rod. The rod is represented as an elongated ellipsoid [20]. In this case, the potential is determined by the expression: components of the electric field strength vector: x, y, z -Cartesian coordinates (i, j, k)-th computational grid node.
The coefficients k x and k z of the nodes surrounding the bar and lying on its axis ((i r , j r , k r ), (i r -1, j r , k r ), (i r , j r , k r -1)), and the coefficient k y of the node at the top of the bar are determined using in (3) instead of U 0 the potential difference at the nodes on the rod and spaced one step of the computational grid from it and using this expression in the formulas for the components of the EF strength (4)- (6). So, for the component E x of the electric field strength, we express x    / through U 0 in the form of the potential difference at a node on the rod axis φ(x ir , y jr , z kr ) = U 0 and at a node spaced from the axis at the distance of the computational grid step φ(x ir-1 , y jr , z kr ), corresponding to the derivative step back: . Let us write the expression for U 0 in terms of considering that the sign of the derivative is taken into account automatically when solving the system of equations: where | ) , For derivatives z    / and y    / |y j =y jr max expressions are similar: The coefficient k xx characterizes the non-linear potential drop between the nodes located on the rod axis and at a step away from it. To determine it, we write Eq. (4) for the cell surface S yz at a distance of half a step from the rod axis (x = x ir+1/2 ) and substitute U 0 into it in the form (7): Then Similar transformations made it possible to obtain expressions for the coefficients k yy , k zz : The coefficients k x , k y , k z (see (2)) are obtained by integrating k xx , k yy , k zz over the corresponding surfaces S yz , S xz , S xy (see fig. 2) using any standard subprogram.
Because the computational domain refers to open areas, then in order to reduce its dimensions when calculating the EF, uniaxially perfectly matched layers (UPML) were placed on its boundaries [17].
The unknown potentials of the nodes of the computational scheme are determined by solving a system of difference equations composed of equations of the form (1) for each node, using the iterative method of alternating directions and the sweep method [16,17]. The boundary conditions for the scalar potential are shown in Fig.1.
The estimation of the accuracy of calculating the EF in systems with rods was carried out by comparing the analytical and numerical solutions for two cases: the potential U 0 is applied to the rod and the rod is in a uniform external EF with the strength E 0 . The analytical solution is made according to the formulas for the potentials and strengths of an elongated ellipsoid [20]. The calculations showed up that for the considered rods with a length-to-radius ratio L/R>10 2 -10 3 , which are under strength, as well as for rods located in an external homogeneous electric field, the maximum values of the relative errors δ Е when calculating the electric field strengths do not exceed 3 %. When calculating based on the approach described in [18], δ Е can reach 45 %, and when calculating based on the approach described in [19], 27 %. Calculations have shown up that a twofold decrease in the step of the computational grid, as well as an increase in its dimensions twofold, do not lead to a change in the distributions of the EF and the values of the relative errors.
An example of electric fields distribution calculation in the vicinity of a rod lightning rod to analyze the conditions for the appearance of upward leaders from it.
The development of positive upward leaders from ground objects occurs when there is an EF strength of the order of 500 kV/m above these objects, caused by an approaching downward leader of negative polarity [4]. The calculation of the EF distribution was carried out for several cases in order to determine the influence on the occurrence and development of upward leaders of a stepped downward leader channel of lightning with a potential of 30 MV-100 MV approaching the object. The boundary conditions used in this calculation are shown in Fig. 1. The case is considered when the space charge is still absent in the zone where the development of the negative streamer stage begins [4]. The step of the spatial grid is chosen equal to 1 m, the absorbing layers located at the boundaries of the computational domain have the number of steps N=10, m=3, k max =300 [18,19].
In the first calculated case, the downward lightning leader is at potential U t and is located at a height of several hundred meters above a grounded lightning rod of length L and radius R. The leader's EF does not affect the EF of the lightning rod top due to the significant distance between them. EF strength in the vicinity of the lightning rod under thunderstorm conditions is E 0 =U t /H t , where H t is the thunderstorm cloud height. The calculation was carried out for the potential of the downward leader U t = 100 MV, launched from a thundercloud from a height of H t = 2000 m. The strength around the top of the lightning rod with a height of L=60 m and a diameter of 2R=0.1 m (see Fig. 1) was E 0 =50 kV/m, the diameter of the lightning leader channel is 2R L =0.01 m. The calculation showed up that the obtained strength above the top of the lightning rod (no more than 100 kV/m) is not enough to fulfill the condition for the appearance of an upward leader. The case was also calculated when the top of the downward leader approached the top of the lightning rod at a distance of 20 m. The distance between the axes of the leader channel and the lightning rod was D LR = 5 m. The strength obtained in this case (300 kV/m) is also insufficient to fulfill the condition of occurrence an upward leader. The calculation under the condition of the location of the downward leader streamer top at a height of 10 m from the top of the lightning rod top showed up the level of the obtained tension exceeding 500 kV/m, thus only in this case the occurrence of an upward leader from the lightning rod is possible.

Conclusions
The proposed method for calculating the EF of conducting rods, in which the non-linear decrease nature of the strength and potential near the rod is taken into account using analytical expressions describing the EF of an elongated conducting ellipsoid under potential, made it possible to reduce the relative error in calculating the EF strength to 3 % or less for rods with L/R> 10 2 -10 3 . In this case, the space step is chosen proportional to the length of the rod, and not to its radius.