Applications of the Parallel-LN-FDTD Method to Calculating Transient EM Field in Complex Power Systems

The present work shows the results of a R&D project carried out by ELETRONORTE and by the Federal University of Para (UFPA). Its core objective is the development of a computational system, called LANE-MAXWELL, for performing analysis and synthesis involving electromagnetic interference (EMI) in a power system substation environment. For the analysis stage, the numerical solution of the problem is obtained by numerically solving the Maxwell’s equations written in a local non-orthogonal coordinate system by employing the Non-orthogonal Curvilinear Finite-Difference Time-Domain Method (LNFDTD). The truncation of the analysis region is done by a new formulation called LN-UPML which involves the solution of the Maxwell Equations for lossy anisotropic media in a general coordinate system. For the synthesis stage, techniques such as neural networks, genetic algorithms and particle swarm optimization are used. The computational environment conceived has a friendly data input/output interface for the user, which permits a better understanding of the electromagnetic phenomena involved. For illustrate the versatility of the computational environment some practical applications involving complex power systems structures are presented.


Introduction
It is widely known that the Finite-Difference Time-Domain (FDTD) method is adequate for solving numerically Maxwell's Equations in order to obtain accurate transient solutions regarding lighting discharges occurring in the vicinity or on power systems (Tanabe, 2001;Oliveira & Sobrinho 2009]. These lighting pulses can produce considerable power in bands 7. Location of PD sources on high voltage (HV) coaxial cables from the obtained transitory signals by using optimization techniques. The necessity of characterization of the electromagnetic phenomena with a greater accuracy led ELETRONORTE to establish an R&D project in cooperation with UFPA (Federal University of Para). This level of accuracy makes the analysis and synthesis of the problems more complex making it unviable by analytic means, in such way the numerical solution techniques became the most adequate choice for the treatment of practical problems. Nowadays, numerical solutions are considered to be as important as experimental measurements, as long as they are complementary to each other. The main objective is to make the practical measurements to have a more economic and safer implementation and procedures. In this context, this work contributes in several ways, based on the development of a computational environment that permits the realization of the analysis and synthesis of problems in electrodynamics in a friendly way. The methodology used in the analyses involves the solution of the Maxwell's equations in Time Domain, written in curvilinear coordinates, by the Finite-Difference Time-Domain method (Taflove & Hagness, 2005). The analysis scenarios are truncated by the uniaxial perfect matching layer technique (UPML): a specific mathematic formalism had to be developed in order to adequate this technique to the curvilinear system (Oliveira & Sobrinho, 2007); for the synthesis several optimization techniques are available, such as: genetic algorithm (Rahmat-Sami & Michielsssen, 1999), particle swarm optimization (Lazinica, 2009) and artificial neural networks (Hagan, & Menhaj, 1994). The software is complemented by input/output interfaces that facilitate the creation of scenarios, visualization of electromagnetic field distribution, video generation, visualization of scenarios, generation of voltage and current graphics, etc.

Related Theory
A. The General Coordinate System Considering a general or curvilinear coordinate system, a position vector r  can be obtained as a function of the general coordinates, and its differential length vector r d  is given by in which the vectors l a  are called unit vectors and form a unit base; which defines the axes of a general curvilinear space (each point in space has its own coordinate system). Figure 1 shows the vectors l a  (l = 1, 2, 3) which are tangent to the l u axes and can be written as functions of the Cartesian Coordinates: x, y and z. An alternative set of three complementary vectors l a  , reciprocal vectors, can be defined in such a way that each one is normal to two unit vectors with different indexes, forming a reciprocal base. This set can be mathematically calculated by the expressions , , , www.intechopen.com in which g is the volume of the hexahedra formed by the vectors , in which g is the determinant of the metric matrix or the covariant tensor (Taflove & Hagness, 2005). Based on the previously presented idea, let a fixed vector F  to be represented by the unit system or by the reciprocal system as shown by (3) .
In (3) . , These components (covariant and contravariant) can be related by means of the following expressions: . , It is worth mentioning that the unit and reciprocal vectors do not necessarily have unitary amplitudes because they depend on the nature of the curvilinear coordinate system used (cells' edge lengths). Hexahedrons similar to that in Fig. 1 in which l F represents the physical components in the base system and their values are thus calculated by (8) . , This representation of vectors can be used to describe electric and magnetic fields' components, as shown in the next topic.

B.
The LN-FDTD Method applied for solving Maxwell's Equations The Maxwell's equations in differential form can be written for lossy, isotropic, nondispersive media as: in which E  is the electric field strength vector, H  is the magnetic field strength vector, ε is the electrical permittivity of the media, μ is the magnetic permeability of the media and σ is the electrical conductivity. By calculating the contravariant components of the fields E  and H  , by using central differences approximations for the derivative terms and by observing the primary and secondary cells in Fig. 2 from which it is easy to see that by observing that, from (2), it is possible to say that � � � . � � � � � �,� , in which � �,� is the Konecker delta function. Equation (11) and the equations for the other five field components can be represented by finite differences by using the standard Yee's algorithm.

C. Stability and Precision of the LN-FDTD Method
The numerical stability of the method is related to the increment Δt, which follows the condition stated by (12), To reduce the effects of numerical dispersion to adequate levels and in order to ensure the precision of the calculations, the dimensions of each non-orthogonal cell in every direction should be less than λ/10 (Taflove & Hagness, 2005), where λ is the smaller wavelength involved in the problem.

D. UPML for Conductive Media -Truncation of the LN-FDTD Method
The Maxwell's equations in frequency domain for an unaxial anisotropic media, can be written as follows: in which  defines the angular frequency of the electromagnetic wave; * E  and * H  are the Fourier transforms of the electric field strength and magnetic field strength, respectively, μ r is the relative magnetic permeability and ε r is the relative electrical permittivity, S is the tensor that promotes the attenuation along the coordinate axes normal to the layers of the absorbing region and σ is the media conductivity. S has the following form www.intechopen.com where S  (α=1, 2, 3) are the parameters that characterize the UPML attenuation along the mentioned general coordinate axes, for which they assume the form (Taflove & Hagness, 2005): By following the procedures similarly to that was indicated in item B (and by introducing auxiliary variables for calculating transformations of fields to time domain), the updating equations for the field components are obtained (Oliveira & Sobrinho, 2007).

E.
Parallel Processing for the LN-FDTD Method The main idea of using parallel computation for solving electromagnetic problems by the LN-FDTD is based in a division of the analysis domain into sub-domains A and B (Fig. 3a). Each sub-domain is a fraction of the whole numerical volume that will be treated by a single processing core. Each core executes essentially the same FDTD code, but with particular boundary conditions. Field information must be exchanged at the sub-domains' interfaces, as illustrated by Fig. 3b. As it is also illustrated by Fig. 3(a), some structures can be divided among two (or possibly more) sub-domains. This issue was treated by implementing a subroutine, which automatically distributes the media parameters among the processing units in such way that the electromagnetic behavior of the divided scenario is identical to the behavior of the full numerical domain processed by a single core. This way, the software can generate complex electromagnetic scenarios by automatically defining the parameters for each CPU by using the data defined by the user in a graphical user interface (GUI), which is graphically displayed by a specific visualization tool developed by using the OpenGL library. This tool, associated to the automatic domain division subroutine, significantly reduces the probability of human errors when constructing and simulating the scenarios. Many figures on this work are direct outputs of the 3D-viewer. The developed software also supports the insertion of thin-wires (Baba et al., 2005) and thinplanes (Maloney & Smith, 1992) (considering the Cartesian coordinate system), and ground ionization due lighting discharges (Ala et al., 2008).

Software Applications
This section, results obtained by using the developed computational environment are presented. Such applications range from the analysis of simple grounding systems to highly complex electromagnetic structures, such as the analysis of the effects of lighting surges on power substations. This way, initially a validation of the implementation is shown, by comparing experimental and numerical results for a grounding structure. Then the importance of using the local coordinate system is evidenced for a grounding structure which geometry is not compatible to the Cartesian system of coordinates. Finally, the effects of lighting currents on three complexes are analyzed for three cases: transmission tower, induced voltages on transmission lines for an urban block of buildings and for the structural part of a power substation.
A. Analysis of a rectangular grounding system. In order to validate the implementation of the parallel numerical algorithm previously described, the transient responses of a rectangular electrode due the application of a artificial lighting current, originally proposed and analyzed by (Tanabe, 2001) were calculated. The geometry is illustrated by Fig. 4, which was reproduced in the computational environment. Basically, this problem consists on a rectangular grounding electrode, with dimensions of 0.5×0.5×3.0 m, buried in a soil with electromagnetic parameters  = 2.28 mS/m,  =  0 and  = 50  0 . Additionally, the problem includes two circuits: one for current injection and other for voltage measurement, which dimensions are defined by Fig. 4. This figure also shows the identification of the sub-domains assigned to each processing core of a Beowulf cluster of PCs (regions numbered from 1 to 10), defined automatically by the automatic domain division subroutine. Fig. 5 shows the numerical results obtained by the simulator developed in this work and the results obtained in (Tanabe, 2001)   According to (Mattos, 2004), the influence of the diameter of the conductors is negligible for the steady grounding resistance R (0.0 Hz), which is given by the equation In (15), L is the rods' length and D is the diameter of the considered circumference. For this case, the expression (15) provides R = 36.224 , which agrees very well to the result obtained by using the LN-FDTD method (Fig. 7) at 2.5 s (R = 36.973). This Figure also presents the results obtained for simulations performed by the classical FDTD algorithm, by approximating the conducting circumference by staircase. It is observed that, as the FDTD spatial step  is reduced, the results tend to converge to the result obtained by the LN-FDTD method (for both, the steady state and transient periods). In Fig. 7, it is also possible to observe that the ohmic values calculated by the FDTD method are higher than the values obtained by the LN-FDTD simulator. This is consistent, as the stair case approximation for the circular path for the electrical current acts as an obstruction for the movement of the electrical charges, thus increasing the voltage/current ratio. This effect is greatly reduced by the non-orthogonal coordinate system and it clearly shows that there is the necessity of avoiding this kind of geometric approximation for performing such analysis especially for higher frequencies.

C. Transmission Line Tower and Induced voltages on Line Insulators due to Remote Lighting Discharge
One of the several applications of the developed computational environment is the simulation of atmospheric discharges and the calculation of induced voltages at specific parts of a structure of interest. For this problem, this structure is a transmission line tower and its line insulators. www.intechopen.com meters away from the tower. The simulations were carried out by considering the following parameters: numerical volume: 250x70x100 m by using cubic Yee cells with edges of 0.5 m; ground parameters:  = 50 0 and  = 2 mS/m; insulator parameters: the insulator modeled has 1.0 m of height with the electrical conductivity of 10 -11 S/m and relative electrical permittivity of 7.5; excitation source: the same used previously (Tanabe, 2001). The lighting discharge was included into the analysis domain by creating a vertical conductor, which penetrates the field absorbing boundary region (UPML) at its top (with impedance matching), representing an infinitely long conductor. By forcing the magnetic field around this conductor to follow the function defined by (Tanabe, 2001), a punctual transient current (52 meters from ground surface) with maximum value of 10 kA was created and used as source. The transmission lines and the grounding rods also penetrate the UPML region. The basis of insulators 1 and 2 (Fig. 8b) were placed at 43.5 meters from the ground surface and the basis of insulator 3 was placed 49.5 meters from the ground plane. It was also simulated the case in which the lighting current occurs at the tower structure, at the point indicated by Fig. 8b and the obtained results for induced voltage on insulator 1 are compared in Fig. 9. Fig. 9 shows the behavior of the voltage induced on the insulator 1 (z-oriented path). In Figure 19, it is shown the induced voltage during the transient period and in Figure 19b it is shown the induced voltage for the stationary period. For this insulator, induced voltage reaches no more than 450 kV during the transient time evolution. As expected, by considering the direct lightning current injection at the tower top, produces a much higher peak, which reaches approximately 1.8 MV at 0.5 s, as shown by Fig. 9a. As shown by Fig.  9b, the difference between the steady induced voltages is approximately 24 kV. For each 10 kA of lighting current peak, the insulator can be subjected to 50 kV (which is added to the nominal line voltage) during considerable amounts of time.

D. Distribution Line in Urban Environment: lighting discharge at a Radio Base Station (RBS).
In recent years, due to world-wide expansion of cell phone systems in urban areas, the raising of tower-like structures operating as radio-base stations (RBSs) became common in www.intechopen.com cities. With considerable heights, around 50 meters, these structures are preferential points for atmospheric discharge strokes. The discharges can occur at the tower's lightning rod or at its metallic body. The occurrences of discharges on telephone towers cause an undesirable set of effects, either in the ground (such as increase of potentials and return currents in grounding systems) or at lines for energy transmission (causing outbreaks due to electromagnetic induction). Such outbreaks are characterized by high-voltage peaks, followed by sudden potential differences between pairs of conductors of the low voltage lines, affecting the consumer units with possible serious damage to electric and electronic devices.
The goal of this application of the LANE-MAXWELL software is to study the induced voltages on low-voltage systems and to estimate how the electromagnetic field, generated by atmospheric discharges on high metallic structures might behave by taking into account the complex structures with different materials around the transmission lines. The two low voltage lines in this scenario can be represented by four metallic cables suspended in the air lined up and separated vertically by a specific distance. On both modeled low voltage lines, two connections between the neutral cable and the ground are established. These low voltage cables were positioned in such a way that their extremities penetrate the UPML of the analyzed region. This set of complex structures modeled in this simulation is represented by a realistic Radio Base Station with its metallic container of equipments and its grounding structure, a set of realistic modeled buildings located next to the RBS, and a realistic road, which was positioned parallel to the x direction of the analysis region as shown by Fig 10. The conception of a scenario was carried out by using the following parameters: 1. Nine eight-floor buildings of 24 m height each, with the following dimensions: 24x12.8x24 m; 2. Analysis region with 168x60x64 m; 3. Yee cells with edges measuring  = 0.2 m ( 840x300x320 cells ); 4. Lateral and frontal separation spaces for Buildings: 6 m and 14 m, respectively; 5. The bars were considered as perfect metal conductors; 6. Wall and street material: reinforced concrete (=0.2 S/m, εr = 7.5 and  r = 1); 7. Low voltage lines with three cables as phase and one as neutral conductors, located at 10 m and 20 m away from the geometric center of a cell phone radio base station (RBS); 8. Cable heights of 6.6 m, 6.8 m, 7.0 m (the three phase lines) and 7.2 m (the neutral line) respectively from the ground. Cables' diameter: 15 mm. 9. RBS is 50.0.m-high and the container has dimensions of 6.0x4.0x3.0 m ; 10. Neutral grounding with ½"x 3 m length cooper rods, with a grounding resistance of 80 Ω. 11. Conductive earth with electrical parameters  = 0.0040836 S/m and ε r = 10. 12. Excitation source: see Fig. 12a. Figure 11 shows the scenario built in the LANE-MAXWELL environment, detailing the above specifications. he obtained trans e RBS (Fig. 12a)  the yz-plane of the structure at x = 84 m (midpoint of the energy line) after 3 μs (12982 iterations) of simulation. It is clear that the electric field which penetrates the buildings is less intense than those present in regions outside this structure. This fact is due to the metallic and concreted constitution of the building structure, which acts as a Faraday's Cage the field, keeping most of the energy of the scattered electromagnetic wave outside the buildings' structure.
(a) (b) (c) Fig. 11. (a) top overview of the urban block, (b) dimensions of one of the buildings and (c) tower, grounding system and transmission lines geometric configuration.
The analysis of the obtained results show high levels of overvoltage and potential differences between live to live and live to neutral conductors on electrical distribution systems, during lightning strokes on a RBS. The results computed with the FDTD method, include full wave solutions for complex structures, involving dielectric materials, a complex tower model, transmission lines, grounding systems, and the consequent complex electromagnetic interactions involving surface waves, diffractions, refractions, and reflections. It has been shown that all those aspects substantially affect results, which, in this work, are close to realistic situations. The problems observed here (relative to transitory high voltage induced in transmission lines) represent high risk when considering the final consumer unit, providing both material and human health damage. Besides that, it was possible to verify from the electric field spatial distribution that high potential differences induced on the soil surface can be found in regions near the tower and its grounding mesh. Such structures are usually located near residences, or in places with constant circulation of people. Telecommunication towers require special attention by the operators which are responsible for the service and by the local energy company as well, keeping safe areas located near those installations, as well as the energy consumer unit.

E. Analysis of Induced voltages on Transmission Lines of a Power Substation due to
Atmospheric Discharge. Figure 14 shows an overview of a computational model of the power substation created by using the LANE-MAXWELL software. It is composed by the structural elements of highvoltage switches, insulators, transmission lines, towers, fence, grit layer, grounding grid, current Transformers, voltage transformers, lightning arrester, etc. The analyzed region's volume is 1000x195x45 m 3 , which is represented by cubic Yee cells with edges measuring 0.5 m. As indicated by Fig. 14, the lighting current was injected at a lighting arrester placed near the 230 kV line, with a maximum value of 1.0 kA (Tanabe, 2001). The wave propagation phenomenon was analyzed up to 12 s after the beginning of the lighting process. The electrical parameters used for representing the ground region were  = 50  0 and  = 2 m S/m. The grit layer was modeled as an isotropic media with thickness of 0.5 m and with  = 0.333 m S/m  = 50  0 . This model follows the specifications of a real substation, which is under responsibility of Eletronorte, installed at the city of Belém (Brazil). In order to create the scenario of Fig.14 with precision, specific database software has been developed. The idea is to create prototypes of the objects of the substation and to clone them at proper spatial coordinates. These objects (towers, switches, etc.) are composed by basic elements of the simulation software (metal and dielectric blocks and thin wires), which are grouped forming a prototype. A text file containing the positions of each clone of each prototype defines the electromagnetic environment. Then, the domain decomposition subroutine distributes the basic elements among the processing sub-domains, which quantity of cores to be used is defined by the user (Fig.3a). Here, sixteen 64-bits processors were used as a Beowulf cluster for solving this problem, which requires approximately 12 hours for concluding the computer simulation. Each core processed a volume of 62.5x195x45 m 3 each time step. Figure 15a shows two y-aligned lines (in doted yellow style) over which induced transient voltages, calculated from the transmission lines to ground, were registered. The red points numbered as 370 and 347 indicate the discrete Yee's indexes related to the x direction, where the voltages are measured. Figure 15b defines of the path used for integrating electric field,    Fig. 16a), it is possible to see that the decay pattern is close to an exponential function with oscillations for the first bay lines. The voltages for the second bay are also close to 10 kV for each kA injected. Figure 16c shows similar results obtained for x = 316. It is possible to see in Fig. 15 (a) that between this point and the discharge coordinates there are several structures, like towers and switches, which act as reflective objects to the electromagnetic field, what justifies the considerable reduction of the peak voltage for the first line (to about 50 kV). The induced voltage to the neighbor line drops to 25 kV and decays almost exponentially to about 6.5 kV at the second bay. The reduction due the reflection of the fields are confirmed when results of Fig. 16d are analyzed. The reference point is x = 1173 (Fig. 15b), which is outside the substation (401 m from the discharge point). The maximum induced voltage is close to 72 kV, decreasing exponentially to 5 kV as distant lines are analyzed. It is also confirmed by analyzing Fig. 17, which shows the electric field distributions for (a) t = 0.225 s, (b) t = 0.365 www.intechopen.com s and (c) t = 1.00 s. In those figures, the blue color represents small amplitudes, green represents intermediate strengths and red is associated to the regions of greater intensity for electric field. For all cases, it is possible to see that the electric field presents smaller intensities as one moves from the fist bay (where the discharge occurs) to the second bay (in y-direction). However, when one moves from the first bay to outside of the substation (in -ydirection), the field attenuation is less intense. Fig. 18 (a) shows the obtained step voltage between the points A and B indicated at Fig. 15b. It is possible to see that after approximately 1 s of the beginning of the lighting discharge, this voltage reaches its maximum value of about 807 V for each kA of peak current. However, its time duration is relatively short. However, the person can be exposed to voltages around 200 V for each kA of lighting current for a much longer time, depending on the time and amplitude characteristics of the injected current. Finally, Fig. 18b shows the electric field distribution at the ground grid plane after 38.5 ms of simulations, illuminating the grid. It is possible to observe that there are red dots, indicating higher intensities of the field, which are associated to the electrical connections of the substation devices to the ground grid. This indicates that current is returning or it is being injected to the grounding mesh.

Conclusion
This R&D project generated a computational environment for the analysis and synthesis of problems on EMC. The initially obtained results were consistent with those available in the literature and with the physics that the problems involve. The used methodology represents a full wave solution and can be used in any frequency range. The developed software can be used in a great range of different applications, as Maxwell's equations are solved by an automated parallel processing environment. Step Potential (V / kA)

Time (s)
Step Voltage (V) www.intechopen.com The association of the computational environment with the laboratory equipment acquired represents a desirable infra-structure for the realization of high level works.