# The Development of Models for Assessment of the Geomagnetically Induced Currents Impact on Electric Power Grids during Geomagnetic Storms.

I. INTRODUCTIONThe functioning of modern electrical power grids is connected with considerable difficulties due to the structural complexity of generating capacities and the basic power network, diversity of their operating regimes, necessity to meet the requirements of reliability and uninterruptedness of electric power supply to consumers, strong external technological and economic ties, uncertainty of future conditions for electrical power grids development, a risk of possible extreme conditions in the system development and other important factors. As a result, modern electrical power grids become more vulnerable to external disturbances, including geomagnetic storms (GMS) [1-5].

During a geomagnetic storm, the variations of the geomagnetic field induce a slowly changing electric field on the Earth's surface. The horizontal component of the geoelectric field is characterized by the intensity of 1-20 V/km and the change period from 10 s to 30 min, the electric-field vector of this field is oriented mainly towards the meridian [6, 7]. EMF appears between the grounded points of power transformers (PT) neutrals of electrical power grid substations [8]. This EMF can reach several kilovolts in long-distance electric power transmission lines (PTL) of 400-500 km, and the quasi direct current, which is also called geomagnetically induced current (GIC), circulates in electric networks [9]. The main impact of GICs on an electrical power grid is the saturation of the power transformers magnetic system [4], [10-13]. While going through the grounded neutrals of the power transformer, GICs produce additional unilateral magnetic biasing of the power transformer core. Consumption of reactive power in the electric power grid increases and transmission capacity of the power network goes down, the voltage in the network drops and false triggering of relay protection and automatic equipment is possible, which can result in accidents and interruption of electricity supply to consumers [4-5].

The studies indicate that vulnerability of electrical power grids during GMS has recently increased due to big geographical areas they occupy and their interconnection and saturation with electronic equipment. The systems of electric power supply and technological processes of consumers are constantly becoming more complicated, economic consequences of the accidents resulting from geomagnetic storms are increasing. Possible interruptions of the electric power supply in major regions can lead to catastrophic consequences for consumers, ecological disasters resulting from the destruction of chemically hazardous domestic and industrial objects [5], [14].

The dynamics of processes occurring in power grids during geomagnetic storms contain multitude of interactive elements and has not been virtually studied. It is practically impossible to carry out a real experiment at an electric power system scale. Therefore, the development of alternative methods and means is needed. One of such alternative methods is modeling that allows to study the regimes of electric power grids (EPG) work under the influence of GICs with values typical of those during geomagnetic storms in the middle and high latitude regions.

II. DEVELOPMENT OF MODELS

The calculations were performed in MATLAB and Simulink environment [15-16]. The developed algorithm of calculation of an electric power grid during geomagnetic storms is shown in Fig. 1. It allows performing the model of electric power grids of any configuration and complexity and those containing both linear and nonlinear elements.

Both Simulink and SimPowerSystems standard blocks and the blocks designed by the authors were used for development of the electric power grid model [17-22]. The initiation of the model takes place every time a calculation of an operation regime of the electric power grid during GMS is made. The initiating involves the analysis of the electric power grid topology, the calculation of the steady state mode and the construction of an equivalent model which then is calculated in Simulink.

A State-Space block or an S-function block is used for modeling of the linear part of the system in the equivalent Simulink model. Simulink models from the powerlib_models library and the models developed by the authors are used to model the nonlinear SimPowerSystems blocks. The nonlinear Simulink models are joined with the linear part of the circuit with the help of input (voltage) and output (currents) variables. The nonlinear models are placed in the feedback circuit of the linear part of the Simulink model. The principal model based on a controllable voltage source or a controllable current source is used for nonlinear elements.

The calculation of an electric power grid model starts with the calculation of a steady state mode in the absence of. In order to do this, prior initiation of the electric power grid is to be performed: the balances of active and reactive powers in the model are to be calculated.

The following special features of the calculated model of the electric power grid, which take into account the GIC flow through its elements during geomagnetic storms, should be noted.

1. The binding of objects to geographical coordinates of a territory is an important stage of the model development as it makes possible to determine the matrix of geoelectric field potentials in nodal points [U.sub.ij] which is used for the calculation of geomagnetically induced currents in power transformer neutrals and windings and in power transmission lines [23].

To record the matrix of geoelectric field potentials at nodal points [U.sub.ij], the differences of potentials on the Earth's surface between each pair of geographical coordinates Ni and [N.sub.j] of grounding conductors of the two substations that are connected by straight portions of overhead transmission power lines are determined with the equation

[mathematical expression not reproducible] (1)

where [mathematical expression not reproducible] is the vector of horizontal component of geoelectric field strength; [mathematical expression not reproducible] is the vector, equal to the length of [dl.sub.j] element of the path and tangential to the chosen path of integration [l.sub.j]; [[alpha].sub.j] is the angle between [mathematical expression not reproducible] and [mathematical expression not reproducible]

Strictly speaking, the formula (1) is valid only if the power transmission lines are located at the same height above the earth surface. It does not take into account the wire sag which has virtually no influence on the result of the calculation due to the fact that the wire sag (up to ten meters) is many times less than the span length between the towers (hundreds of meters) and the overall length of PTL (tens, hundreds of kilometers).

It should be noted that the geoelectric field [mathematical expression not reproducible] is the rotational one. However, the path of integration in the formula (1) can always be superposed with the route of PTL. If we suppose the geoelectric field to be practically uniform (E=const) within the limited geographic area occupied by EPG and if we divide the route of PTL into rectilinear sections, it can be possible to substitute the integration in the formula by the summation

[mathematical expression not reproducible] (2)

where [l.sub.j] is the length of the [j.sup.th] sector of the power transmission line; [[alpha].sub.j] - orientation angle of the [j.sup.th] sector of the power transmission line relative to the electric vector of the geoelectric field; [L.sub.E] is the electrical distance between the substations that is equal to the sum of the projections of the power transmission line rectilinear sections on the geoelectric field vector direction.

Despite the given limitations, the formula (2) can also be generalized in the case of the non-uniform geomagnetic field. In order to do so it is enough to divide the geographic area occupied by the EPG under consideration into several similar segments, for instance, of rectangular shape, within which the electric field is uniform. Then the formula (2) will remain true within each segment.

2. GIC matrix for the power network consists of GIC values for every [ij.sup.th] closed circuit of their flow (tellurics--a grounding lead--a neutral--a power transformer grounded winding [T.sub.i]--a PTL wire - a power transformer grounded winding [T.sub.j]--a neutral - a grounding conductor--tellurics) and is written like:

[mathematical expression not reproducible] (3)

where [mathematical expression not reproducible] is the total resistance matrix for the GIC flow

contour. In order to calculate separate elements [I.sub.GICij] of the GIC matrix with the help of the formula (3) one should use the corresponding element [U.sub.ij] of the matrix of potentials which has been preliminary estimated by the formula (2) and the element [R.sub.[SIGMA]ij] of the total resistance matrix which is estimated with the help of the expressions

[mathematical expression not reproducible] (4)

where [r.sub.Lij] is active resistance of the phase wire of PTL; [r.sub.HVj], [r.sub.HVi] - active resistance of the high-voltage winding of the power transformers [T.sub.i] and [T.sub.j] respectively; [r.sub.Gi], [r.sub.Gj] - resistance of grounding grids of the power transformers [T.sub.i] and [T.sub.j] respectively.

3. To calculate the magnetizing currents the power transformers of electrical power grids are modeled on the basis of application of a T-shaped equivalent circuit with the nonlinear inductance of the magnetizing branch. In order to calculate the magnetizing currents in the power transformers the following conditions have to be set

- the structural parameters of power transformers (the section of magnetic system rod, average length of magnetic line, number of turns of the grounded windings);

- the parameters of the magnetic hysteresis loop of the power transformer system steel: the coercive force [H.sub.c], the residual flux density [B.sub.r] and the saturation flux density [B.sub.s] for the magnetization curve modeling;

- the parameters of the equivalent circuit of the electrical power grids;

- the matrix of the geoelectric field potentials at the nodal points of the electrical power grids.

Magnetization curve of the power transformer magnetic system in the model is set as the piecewise-linear dependence between flux linkage [psi] and instantaneous magnetizing current value [i.sub.0]. The dependence [PSI]([i.sub.0]) is determined from the B=f(H) curve model with the help of a simplified hysteresis model of the power transformer magnetic system [25]. In order to build [PSI]([i.sub.0]) curve 150-200 values are used in the calculated model of the power transformer.

When geomagnetically induced currents flow in the grounded windings of the power transformers, the inductance of the magnetizing branch L=[PSI]([i.sub.0]) / [i.sub.0]=L([I.sub.GIC]) is a function of the GIC value and it is determined by the location of the working point on the magnetizing branch.

Instantaneous values of the magnetizing current of power transformer [T.sub.i] that is a part of the [ij.sup.th] closed contour of GIC flow are determined from the solution of the following nonlinear differential equation

[mathematical expression not reproducible] (5)

where [u.sub.i](t) is the voltage of the magnetizing branch of the power transformer [T.sub.i] winding; [i.sub.0i]. is the magnetizing current of the power transformer [T.sub.i]; [L.sub.i]([I.sub.GIC]) - inductance of the magnetizing branch of the power transformer [T.sub.i].

4. When GICs flow in grounded windings of the power transformer, the peak values of the magnetizing currents increase repeatedly and have a notable half-wave asymmetry. The power transformers magnetizing current becomes the source of the high-order harmonics in electric power grids during geomagnetic storms. Therefore, to calculate the instantaneous values of the current and voltage in the elements of the estimated model of the electric power grid, the harmonic component of the magnetizing currents [i.sub.0i](t) of the power transformer is preliminary determined [mathematical expression not reproducible] (6)

where [I.sub.0i(n)], [[phi].sub.0i(n)] are the amplitude and the initial phase of the magnetizing current n-th harmonic of the power transformer [T.sub.i]; and [omega] is the angular frequency of the system voltage

5. The processes of the power transformer magnetic system saturation are accompanied by multiple increases of the non-sinusoidal magnetizing currents that propagate in the electric network which results in a considerable distortion of the currents and voltages curves in the power grid. Thus, active power of the non-sinusoidal current in an element of the calculated model of the power grid is determined as an average power for the period [mathematical expression not reproducible] (7)

[mathematical expression not reproducible] (8)

where [U.sub.(n)], [I.sub.(n)] are the RMS values of n-harmonics of voltage and current; [[phi].sub.(n)] = [[phi].sub.u(n)] - [[phi].sub.i(n)] The apparent power S is defined as the product of RMS currents and voltages:

[mathematical expression not reproducible] (9)

The reactive power in an element of the calculated model of the power grid is defined as the sum of reactive powers of separate harmonics:

[mathematical expression not reproducible] (10)

Thus, the following steps were realized in the algorithm of calculation of the electric power grid regimes during electromagnetic storms:

- binding of objects of an electric power grid to the geographical coordinates of the territory;

- calculation of geomagnetically induced currents in neutrals and grounded windings of the power transformers, in the phases of PTL;

- calculation of the instantaneous and RMS magnetizing currents of power transformers of the power grid;

- decomposition of the currents sources of high-order harmonics for modeling of the power transformers magnetization branches into constant and sinusoidal components, i.e. Euler-Fourier series;

- calculation of instantaneous and RMS currents and voltages for every harmonic in EPG elements;

- calculation of the RMS active, reactive, apparent capacities with the consideration of the impact of harmonic components of the current and voltage in EPG elements.

III. SIMULATION RESULTS

In order to assess the GIC impact on an electric power grid a model of electric power grid of Samara region was developed. Characteristic features of Samara region EPG are inherent to any high-power EPG: availability of major power sources--thermal power station "TPS of AVTOVAZ", "Togliatti TPS" and "Zhiguly Hydroelectric Power Plant" with total generating capacity of 4449 MW; availability of the intersystem and intrasystem communication lines with voltage 110, 220 and 500 kV and total length of over 1000 km; availability of the major nodal transformer substations with voltage levels of 500/220/110, 220/110/10 kV. In the model of EPG the substations for power supply of industrial and urban load with voltage levels 220/10, 110/10, 110/6.3 kV are taken into account. The total installed capacity of the power transformers and autotransformers at the power stations and the transformer substations in the EPG model is 10883 MVA. The structural network of Samara region EPG model is shown in Fig. 3.

The power transformers and autotransformers of EPG are implemented with grounded neutrals. There is a big number of long EPLs of all voltage in the EPG. The EPLs are mainly west-east oriented. Therefore, during an intensive geomagnetic storm and when the direction of the geoelectric field coincides with the direction of EPG objects location on the geographical map of the territory, considerable geomagnetically induced currents are possible in the neutrals and grounded windings of PT and phase wires of PTL.

In the calculated model:

- the pulse shape of geoelectric field intensity during a geomagnetic storm is shown in Fig.4 and corresponds to the pulse shape in the work [26], the maximum values of intensities are assumed equal to 6, 10, 15, 20 V/km;

- the west-east direction of the geoelectric field strength vector E was chosen;

- the calculation of regimes of the EPG work during geomagnetic storms was made from 300 sec to 1300 sec. The impact of the geoelectric field intensity vector on the power grid starts in 370 sec. from the beginning of the model calculation. The peak value of pulse is at the 1280-th sec.

The results of computer simulation for various geoelectric fields are shown in Fig. 5 - 8.

The shape of the curve of geomagnetically induced current in HV winding of the substation power transformer repeats the shape of geoelectric field impulse curve, the maximum value of GIC is defined by the maximum value of the strength of geoelectric field E and varies from 15.4 A for E = 6 V/km to 82.6 A for E = 20 V/km (Fig. 5).

The peak values of the magnetizing current several times exceed the values of the no-load current of the power transformer: when E = 6 V/km the peak values of the magnetizing current reach the value of 151 A, when E = 20 V/km - 645 A. The magnetizing current curve is strongly distorted, the effect of half-wave asymmetry is observed (Fig. 6).

The total current curves of HV winding of the power transformer are also strongly distorted, there are the 2nd, 4th and constant harmonics of the current (Fig. 7).

The line spectra of harmonic components of HV winding current of power transformer 110/6.3/6.3 of "Elkhovka" substation for different intensities of geoelectric field are shown in Fig. 8.

It should be noted that the nature of the current harmonics in HV winding of a power transformer is stable independently of the geoelectric field intensity. The DC component is 7.5-21.5 % from the fundamental harmonic (E = 6-20 V/km). The total harmonic distortion of the current in HV winding of the power generator is [K.sub.I(n)] = 18.64--37.3 % (E = 6--20 V/km). The total harmonic distortion of the current [K.sub.I(n)] with regard to current harmonics up to the 40TH order in the point of electric power transmission is calculated by the formula:

[mathematical expression not reproducible] (11)

where [I.sub.(1)]--is the current of the fundamental harmonic. The results obtained by the computer simulation correlate qualitatively well with the statistical data on GIC impact on the electric power grids, [4, 6, 27, 28].

IV. CONCLUSION

The developed calculation methods allow defining geomagnetically induced currents in power grids elements with regard to their geographical location on the territory map. The computer-aided implementation of the developed methods in Simulink, an add-on product to MATLAB, with the use of the blocks developed by the authors makes it possible to model processes of power transformer saturation during geoelectric storms and to define the instantaneous and RMS values of magnetizing currents of power transformers; RMS values of geomagnetically induced currents in neutrals and grounded windings of PT and phase wires of EPL, the instantaneous values of total currents and voltages in all the elements of the electric power grid. It also enables to perform the harmonic analysis to assess the geomagnetically induced currents impact on the power grid working regimes, which makes possible to define acceptable GICs values for the reduction of their negative impact.

REFERENCES

[1] "Report of the Commission to Assess the Threat to the United States from Electromagnetic Pulse (EMP) Attack. Critical National Infrastructures", 2008. [Online]. Available: http://www.empcommission.org/docs/A2473-EMP_Commission-7MB.pdf

[2] D. H. Boteler, "Assessment of geomagnetic hazard to power systems in Canada," Natural Hazards, vol.23, Is.2, pp. 101--120, 2001. [Online]. Available: http://dx.doi.org/10.1023/A:1011194414259

[3] J. Watermann, "The magnetic environment--GIC and other ground effect," Space weather, vol. 344, pp. 269--275, 2007. [Online]. Available: http://dx.doi.org/10.1007/1-4020-5446-7_23

[4] J. Kappenman, "Geomagnetic Storms and Their Impact on the U.S. Power Grid," Metatech, Meta-R-319, 197p, 2010. [Online]. Available: http://web.ornl.gov/sci/ees/etsd/pes/pubs/ferc_Meta-R-319.pdf

[5] J. Kappenman, "A Perfect Storm of Planetary Proportions" IEEE Spectrum, vol.49, Issue 2, pp. 26--31, 2012. [Online]. Available: http://dx.doi.org/10.1109/MSPEC.2012.6139230

[6] A. J. McKay, "Geoelectric Fields and Geomagnetically Induced Currents in the United Kingdom" University of Edinburgh, p. 260 2004. [Online]. Available: http://hdl.handle.net/1842/639

[7] A. Pulkkinen, E. Bernabeu, J. Eichner, C. Beggan and A. W. P. Thomson, "Generation of 100-year geomagnetically induced current scenarios" Space Weather vol 10, issue 4, 2012. [Online]. Available: http://onlinelibrary.wiley.com/doi/10.1029/2011SW000750

[8] A. Pulkkinen, "Geomagnetic flux density during highly disturbed space weather conditions: studies of ground effects," Academic Dissertation. University of Helsinki, Faculty of Science, Department of Physical Sciences, Helsinki, 2003 [Online]. Available: http://ethesis.helsinki.fi/julkaisut/mat/fysik/vk/pulkkinen/geomagne.pdf

[9] R. Pirjola, "Review on the calculation of surface electric and magnetic fields and of geomagnetically induced currents in ground-based technological systems," Surveys in geophysics, vol.23, Issue 1, pp. 71-90, 2002. [Online]. Available: http://dx.doi.org/10.1023/A:1014816009303

[10] R. Girgis, K. Vedante and K. Gramm, "Effects of GIC on power transformers and power systems" Transmission and Distribution Conference and Exposition (T&D), 2012 IEEE PES, pp. 1-8, 2012. [Online]. Available: http://dx.doi.org/10.1109/TDC.2012.6281595

[11] P. R. Barnes, B. W. McConnell and J. W. Van Dyke, "Electromagnetic pulse research on electric power systems: Program Summary and recommendations," 1993. [Online]. Available: http://dx.doi.org/10.2172/6721434

[12] R. A. Walling and A. N. Khan, "Characteristics of transformer exciting--current during geomagnetic disturbances," IEEE Transactions on Power Delivery, vol. 6, Issue 4, pp. 1707--1714, 1991. [Online]. Available: http://dx.doi.org/10.1109/61.97710

[13] R. Pirjola, "Geomagnetically induced currents in the Finnish 400 kV power transmission system." Physics of the Earth and Planetary Interiors, vol. 53, Issue 3-4. pp. 214--220, 1989. [Online]. Available: http://dx.doi.org/10.1016/0031-9201(89)90005-8

[14] J. G. Kappenman, "Space weather and vulnerability of electric power grids," Effect of space weather on technology infrastructure, pp. 257 - 286. 2004. [Online]. Available:http://dx.doi.org/10.1007/1-4020-2754-0_14

[15] K. E. Lonngren and S.V. Savov, "Fundamentals of Electromagnetics with MATLAB," Scitech, 2005, 574 p.

[16] "SimPowerSystems User's Guide (Second Generation)", 2013. [Online]. Available: http://www.mathworks.com/help/pdf_doc/physmod/powersys/powersys.pdf

[17] V. V. Vakhnina, V. D. Selemir, V. I. Karelin, V. A. Shapovalov and V. V. Gorokhov "Model of Power Supply System of the City", 2008, March 25, Certificate of state registration of computer programs, RU No. 2008611506.

[18] V. V. Vakhnina, A. N. Chernenko, M. S. Makeev and V. A. Kuznetsov, "Model of a Power line", 2013, February 6.Certificate on state registration of computer program, RU No. 2013611832.

[19] V. V. Vakhnina, A. N. Chernenko, M. S. Makeev, and D.A. Kretov, "Model of Power Transformer Magnetizing Branch" 2013, February 6, Certificate of state registration of computer program, RU No. 2013611833.

[20] V. V. Vakhnina, A. N. Chernenko, D. A. Kretov and V. A. Kuznetsov, "Model taking into account the power transformer magnetizing branch nonlinearity", 2013, April 30, Certificate of state registration of computer program, RU No. 2013615269.

[21] V. V. Vakhnina, A. N. Chernenko and M. S. Makeev, "Dynamic Load Model for Power Supply System", 2012, November 22, Certificate on state registration of computer program, RU No. 2012660528.

[22] V. V. Vakhnina, A. N. Chernenko, M.S. Makeev, and V. A. Shapovalov, "Model of High-voltage Circuit Breaker", 2013, February 6, Certificate on state registration of computer program, RU No. 2013611831.

[23] V. V. Vakhnina, D. A. Kretov and V. A Kuznetsov, "Calculation of geo induced currents in high"voltage power lines of electricity supply systems at geomagnetic storms," Proceedings of the Samara Scientific Center, vol. 14. No. 6.2012. pp. 244--246, 2012. [Online]. Available: http://www.ssc.smr.ru/media/journals/izvestia/2012/2012_6_244_246.pdf

[24] V. V. Vakhnina, A. A. Kuvshinov and V. A. Kuznetsov "Reduction of Risks of Accidents in Power Grids during Geomagnetic Storms" Heliophysical Study, Issue 5. pp.115--123, 2013. [Online]. Available:http://vestnik.geospace.ru/index.php?id=191.pdf

[25] J. H. Chan, A. Vladimirescu, X.-C. Gao, P. Liebmann and J. Valainis, "Nonlinear Transformer Model for Circuit Simulation," IEER Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol.10, Issue 4, pp. 476--482, 1991. [Online]. Available:http://dx.doi.org/10.1109/43.75630

[26] P. R. Barnes, D.T. Rizy, B.W. McConnel, F.M.Tesche, E.R Taylor. (1991); OAK RIDGE NATIONAL LAB TN. Electric utility experience with geomagnetic disturbances. Under interagency agreement No. 0046-[C.sub.1]56-A1, 98 p. http://web.ornl.gov/~webworks/cpr/v823/rpt/51089.pdf

[27] J. Elovaara, "Finnish experience with grid effect of gic's," Space Weather. Astrophysics and Space Science Library, vol. 344, pp. 311 - 326, 2007. [Online]. Available:http://dx.doi.org/10.1007/1-4020-5446-7_27

[28] A. Pulkkinen, A. Viljanen, R. Pirjola "Large geomagnetically induced currents in the Finnish high-voltage power system," Finnish Meteorological Institute, Helsinki, 2000, 99 p.

Vera V. VAKHNINA, Aleksey A. KUVSHINOV, Vladimir A. SHAPOVALOV, Aleksey N. CHERNENKO, Dmitry A. KRETOV

Togliatti State University, 445667, Russia V.Vahnina@tltsu.ru

Digital Object Identifier 10.4316/AECE.2015.01007