Power Quality Blog

Distributed Photovoltaic Integration as Complementary Energy: Consideration of Solutions for Power Loss and Load Demand Growth Problems

Published by Salama MANJANG¹, Yuli Asmi RAHMAN²,
University of Hasanuddin (1), University of Tadulako (2), Indonesia

Abstract. The purpose of this study is to optimize the location and capacity of PV in the feeder distribution system 20 kV of Central Sulawesi, Indonesia. The proposed method uses the optimization method of development from the genetic algorithm, namely NSGA-II. Optimization is carried out in three scenarios by considering the value of the total active PV power capacity which produces the minimum active power loss and voltage deviation. The simulation result shows that the integration of PV-DG can improve drop voltage of distribution system performance due to load growth effect.

Streszczenie. Celem tego badania jest optymalizacja lokalizacji i wydajności PV w systemie dystrybucji zasilania 20 kV w środkowym Sulawesi w Indonezji. Proponowana metoda wykorzystuje optymalizację opartą na algorytmie genetycznym, mianowicie NSGA-II. Optymalizację przeprowadza się w trzech scenariuszach, biorąc pod uwagę wartość całkowitej mocy czynnej PV, która powoduje minimalne straty mocy czynnej i odchylenie napięcia. Wynik symulacji pokazuje, że integracja PV-DG może poprawić wydajność systemu dystrybucji ze względu na efekt wzrostu obciążenia. (Rozproszona integracja fotowoltaiczna jako energia uzupełniająca: rozwiązań problemów związanych z utratą mocy i wzrostem zapotrzebowania na obciążenie).

Keywords: PV on grid, load growth, genetic algorithm, distributed generation, drop voltage
Słowa kluczowe: PV na siatce, wzrost obciążenia, algorytm genetyczny, generacja rozproszona, spadek napięcia.

Introduction

The estimated increase in global energy demand in 2040 reaches 37% [1]. This was triggered by an increase in electricity consumption from the residential and commercial sectors each year. This certainly has an impact on increasing the load on the distribution network and impacting the operational problems of the system. To anticipate this, the use of distributed generator in the distribution system is one alternative solution to overcome the problem of increasing the load and effects caused on the operational side of the distribution network. Distributed generation (DG) is the opposite of the traditional model of centralized generation. If in the centralized generation model a large-scale power plant provides energy for one region, the DG prioritizes many small capacity power plants that are connected or not connected to the distribution and transmission network. In general, DG refers to power plants with a capacity of 100 kW untill 50 MW [2]. Based on International Energy Agency (IEA) statistical data, the average energy loss during distribution and transmission in a centralized electricity generation system is the range between 8 and 15% [3]. This is a challenge that arises from centralized electricity generation but is also an opportunity for DG innovation as decentralized generation; which is closer to the consumer and utilizes local renewable energy.

The trends in global energy use moving towards renewable energy use, including in the electricity generation sector, also affect the transition of various sectors from fossil fuels to electricity utilization. On-grid photovoltaic (PV on-grid) is currently the latest trend in the electricity business model.

In general, Indonesia has a solar energy potential of 4.8 Kwh / m^² , equivalent to 112.999 GW peak (GWP). Palu as the capital of Central Sulawesi province with an astronomical location of 0^o57 ‘LS; 120^o0 ‘BT 5.512 so that it is right on the Equator line with an altitude of 0-700 making this region very potential for PV placement as mapped in Fig. 1. Referring to National Aeronautics and Space Administration (NASA) data, the potential of solar radiation in the city of Palu in 2017 ranged from 4.98 – 5.86 / m² per month with the duration of sunlight 6-8 hours per day shown in Fig. 2.

Fig.2. Solar radiation per month in Palu city

Related research on the potential of solar energy in the city of Palu revealed a comparison of the efficiency levels of the three types of solar cell module technology. The output power of monocrystalline, polycrystalline, and amorphous silicon each varies between 20-23 kWh / m², 16-19 kWh / m², and 7–9 kWh / m². Thus, there is a promising potential for the supply of electrical energy from solar panels as a DG source that inject active power. However, the on-grid PV business is a challenge for providing adequate and appropriate reserve margins to offset the intermittency of PV [5] and the use of environmentally friendly technologies. PV on the grid has a role in the quality of power in the distribution system. The role is to improve the voltage value, improve the impact of load imbalance, reduce the occurrence of active and reactive power fluctuations, increase the power factor, reduce the number of power losses and improve the value of system reliability [6-9]. PV integration in the distribution feeder must meet the requirements of the connection technique, namely the selection of the location (placement) and capacity (sizing) of PV power. Determination of on-grid location and PV power capacity is proposed using optimization techniques. Previous studies have used optimization methods in determining the connection point of DG with the distribution network [10-13].

However, the addition of PV as DG which is generally placed near the load on the distribution network will greatly affect the overall system power flow. A significant influence and needs to be considered further is system loadability. Loadability provides the maximum value for loading that is still possible without disrupting the stability of the system work. The maximum loadability parameter in a distribution system that must be considered is the bus voltage limit. This voltage limit will be a limitation of the optimization that will be done in determining the size and location of DG placement to be installed on the system [14]. To optimize the function of PV as an additional plant in meeting the increasing electricity demand, planning for the selection of the placement location and determination of the DG unit’s output power must be carried out appropriately. Therefore, this optimization problem will be solved by using the NonDominated Sorting Genetic Algorithm-II (NSGA-II) method. This method can solve multiobjective optimization problems in the placement and determination of PV-DG output power so that the formation of the number, location, and PV output power in the maximum loadability state of the system can still be achieved.

Power flow after integrated PV-DG

This research aims to optimize the placement of PV-DG in the system to produce minimal power loss. The calculation to determine the total power loss in the system uses the power flow equation described in the following section.

Fig.3. A Single – line diagram of distribution feeder

Figure 3 represents the distribution system with a current value of:

where

The current value J of equation (1) is obtained by forming the Bus current Injection to Branch Current (BIBC) matrix, that is,

By using Kirchhoff’s voltage law an equation is obtained,

Power loss on each line connecting bus k to bus k + 1 can be stated as follows:

From equation (5), the total power loss can be calculated by summing the power loss in each line so that the following equation is obtained

Non-Dominated Shorting Genetic Algorithm-II (NSGA-II)

Solution optimization problems that involve multiple objectives (multiobjective optimization) have an impact on increasing the number of optimal solutions to these multiple problems, this is widely known as Pareto-optimal solutions. Some evolutionary algorithms for solving multi-purpose problems, Multiobjective Evolutionary Algorithms (MOEAs) have been discovered [15-17].

The main reason for using solutions with this method is their ability to find many Pareto-optimal solutions once a simulation is run. Since evolutionary algorithms (EAs) work with a population that contains problem-solving, a simple EA can be developed to maintain several different sets of expected solutions. With an emphasis on moving towards the actual Pareto-optimal region, an EA can be used to find several Pareto-optimal solutions with a single simulation. NSGA is the first evolutionary algorithm to adopt it.

However, in its use, this algorithm still has many weaknesses in the computational process, the absence of elitism that prevents the loss of a good solution when the solution has been found requires a specific definition of the value of sharing parameters (σshare). Therefore it is necessary to change the original code from NSGA. These changes are referred to as the second version of the NSGA or abbreviated as NSGA-II.

From various simulation results on several multiobjective problems that are difficult to solve, NSGA-II shows very satisfying results [18] in terms of the discovery operation of a group of problem-solving with high variation, as well as the proximity of the convergence with the group of Pareto-optimal solutions.

However, the most important advantage in the selection of NSGA-II is its ability to overcome the three weaknesses of NSGA above, because NSGA-II has a better sorting algorithm, the incorporation of the concept of elitism in the main algorithm, and no need to share parameters that must be determined first. Some stages in the NSGA-II in the outline are; population initialization, non-domination sort, crowding distance, crowded-comparison-operator, selection, genetic operator 1 (crossover), genetic operator 2 (mutation), and recombination [18]. The implementation of the NSGA-II in determining the optimal size and location of the PV is shown in the flow chart of Fig. 4. Fitness calculation for each individual is based on the objective function (7) to be achieved in the placement and determination of the optimal capacity of this PV. The objective function consists of two functions to minimize system power loss and voltage deviation.

Fig.4. Flowchart for PV-DG placement using NSGA-II

As explained in the background of the problem, the placement of DGs is done to improve the system load supply capability, so that it can accommodate the additional load that will occur. Addition to load affects the stability of the system voltage.

Subjects to (7-10) :

Each individual will produce different fitness values. Individuals with each fitness value will be combined in one array and sorted according to the rules of non-domination sorting and crowding distance so that individuals are ranked the best to the worst for the selection process. Some parameters used to generate the initial population are tabulated in Table 1. Each PV unit is represented by S binary strings of 8 gene bits in one chromosome. The first four bits represent the PV unit capacity, while the remaining 4 bits represent the PV output power. Following the PV output power rating, the output power that can be produced is between 1-5 MW. Whereas the DG placement location is encoded by a binary L string of 4 bits. The formation of chromosomes is done by combining the S and L strings into one array of individual genes consisting of 12 bits in which the value of each chromosome element (gene) is generated randomly using MATLAB 2016a software.

Table 1. The parameters of NSGA-II

Result and discussion

Scenario 1 : existing (base condition )

Initially, the Tipo system was analyzed for load flow that occurred in each of these systems without the addition of DG. Power flow analysis is performed using the Newton Rapson method, with accuracy reaching 0.00001, acceleration 1.6, and maximum iteration of 100 times the process.

System parameters to be considered in the simulation are the total network losses and the voltage value of each bus in the system.

Fig.6. The power flow of the existing 13 bus GI Tipo system.

Power flow information on the Tipo system was analyzed with a Matlab based program. The results of the power flow study are shown in Fig. 6 and Fig. 7. The total system load of 19,980 MW + j 9,678 MVAr is served by the system of 21,221 MW + j 10.985 MVAr. Figure 7 shows the value of power flow and power loss in the line. The total power loss that occurred was 1.271 MW + 1.135 MVAr. The biggest power losses are in the Luar Kota line, Kelor line, and Dalam Kota line respectively 0.616 MW + j0.541 MVAr; 0,107 MW + j0,252 MVAr; and 0.294 MW + j0.251 MVAr. Figure 7 shows the minimum voltage conditions on buses connected to Tipo’s Electrical Substation of 0.806 pu experienced by Luar kota bus. The value of the voltage drop that occurs is 16.8%.

Scenario 2: placement of PV-DG

Power flow studies conducted on Tipo system were discussed in the previous section. The results of the study show that three buses have very poor under voltage conditions. Radial topology condition does not allow DG placement of less than three units.

Fig.7. Comparise Voltage bus profil and load bus

Optimization of on grid PV placement by one unit and two units has been carried out. Determine a portion of the amount on the bus. Placement of PV-DG by one unit up 1.2 MW on the Luar Kota bus can reduce active power loss by 5.4% and reactive power loss by 12.73%. Voltage profile of Dalam Kota and Luar Kota buses increased respectively 12.15% and 18.26%. However, the voltage condition on the Kelor bus in conditions under-voltage reaches 0.862 pu.

The same thing happened when the placement of two PV-DG units on Luar Kota bus and Kelor bus of 2 MW and 3.1 MW, respectively, resulted in a reduction in active and reactive power losses of 28.32% and 11.61%. The voltage on the Luar Kota bus and Kelor bus corrected respectively 18.26% and 15.79% but the voltage on the Dalam kota bus is still in an under voltage condition at 0.914 pu. The results of optimization using NSGA-II are shown in Table 2.

Table 2. The results of the optimization of the PV-DG placement on the 13 Tipo’s Electrical Substation bus system

Table 3. Recapitulation of PV-DG placement effect on 13 Tipo’s Electrical Substation bus systems

The integration of 5.7 MW from DG-1 can reduce active and reactive power losses by 61.81% and 48.11% respectively so that the system power loss in the Tipo GI became 0.569 MW + j0.877 MVAr as tabulated in Table 3. The reduction in power losses that occur can be explained as a result of injection of the current generated from the PVDG so that it can reduce the current from the main generator which is centralized.

Scenario 3

The type load model on the 20 kV distribution network Palu is mixed. The industrial load is represented by industrial bus and express bus 6 represents commercial loads. The other buses are buses with resident loads. This scenario begins by testing the DG placement solution for each type until the third growth year.

The load growth rate refers to the national rate of 8.6%/year. In this test, the assumption does not change the distribution network in terms of conductor size or network configuration.

Fig.8. Power loss of Tipo 13 bus GI due to the ‘g’ factor

The test results show that the growth of the load in the first year resulted in 25.794 MW + j12.957 MVAr and followed by an increase in a power loss of 16.65% or to 1.525 MW + j1.689 MVAr. The same thing happened to the load growth in the following year which resulted in an increase in load to 26.456 MW + j12.757 MVAr and 28.207MW + j13.278 MVAr, respectively.

Fig.9. System bus 13 Tipo’s Electrical Substation bus voltage profile due to the ‘g’ factor

The load growth also affect the voltage profile of each bus in the Tipo GI system. In the first year, the biggest voltage drop occurred on the Kelor bus as the farthest bus from the slack bus valued at 16.3%. Such is the case in the second and third years of 15.8% and 17.7%, respectively. The following test is find out the effect of the system up to the third growth year if the PV-DG integrated system is on three predetermined buses. A comparison of power losses and voltage profiles between system conditions without PV and PV integrated systems is shown in Fig. 8 and Fig. 9.

Fig. 8 shows that the integration of PV-DG can minimize power losses that occur until the third year. The first year of the DG-1 integration of 5.7 MW was able to reduce power losses by 6.94% to 1,426 MW + 1,692 MVAr. Likewise in the second and third years, there was a decrease in power losses to 0.807MW + j1.21MVAr and 0.962MW + j1.415MVAr respectively. PV on-grid can increase the minimum voltage value when system conditions are without DG. However, this improvement was unable to maintain the voltage level within the tolerance range of 0.95 pu to 1.05 pu as shown in Fig. 9.

The voltage profile of Dalam Kota bus in the first, second and third years show values of 0.948 pu, 0.936 pu and 0.918 pu. Flamboyan bus also experienced a voltage drop in the third year exceeding the limits of 0.94 pu.

Based on the above analysis, the DG placement scenario on the Tipo system is carried out by taking into account the ‘g’ growth factor. In the optimization process, there are changes in the parameters used in the optimization process, that’s :

• Total PV-DG = four units
• Limit of obstacle : 0,1 ≤ P_DG ≤ 8,5 MW

The change in the amount of SG is based on a test scenario conducted previously using three units whose results are not converging. The results obtained violate the voltage limit. For that scenario, the number of DGs is increased to four units according to the previous optimization recommendations. The results of these scenarios are shown in Table 4.

Table 4. The results of optimization of DG-1 placement on the Tipo system consider the faktor g ’factor

DG’s optimal solutions shown in Table 4 have the effect of reducing power loss in the first, second and third year. The test results show that in the first year load growth with a total load of 25.794 MW + j12,957 MVAr, DG 8.2 MW integration can reduce power losses significantly by 43.14%.

Fig.10. Power loss before and after DG integration (optimization results Table 4) in each growth year.

In the second year, the system load of Tipo system became 26,456 MW + j12,757 MVAr followed by an increase in power loss to 1,758MW + j 1,942MVAr. PV-DG integration can reduce real power losses by 49.94%. In the third year, the Tipo system bears a total load of 28.207MW + j13.278 MVAr. The integration of DG 8.2 MW can reduce power loss by 53.4% of system power loss without DG by 2.077MW + j2.28MVAr. The impact of PV-DG integration on reducing power losses with a three-year scenario of load growth is presented in Fig. 10. In addition to reducing power losses, DG integration is also able to improve the voltage profile that experiences voltage drops on certain bus due to load increases.

In the first year, the biggest voltage drop occurred on the Kelor bus as the farthest bus from the slack bus valued at 14.1%. Likewise in the second and third years respectively 15.8% and 17.7%. PV-DG integration can reduce the value of voltage drop in the first, second, and third years by 1.4% each; 2.6%; and 3.9%.

To test the performance of NSGA-II, simulations with the same procedures and objects have been carried out using Simple genetic Algorithm (SGA). The results show improved convergence generation and computational time trimming with more significant results. This is shown in Fig. 11.

Fig.11. SGA-NSGA-II performance parameters

The results of this study are also supported in previous studies [19] that have tested the placement of PVDG in the IEEE 30 bus test system. There are different methods and objective functions in this study which only focuses on reducing the active power loss of the line. The method proposed in this study is the Breeder Genetic Algorithm (BGA). Similar research has also been carried out [20] with the same method but with a different type of DG.

Conclusion

This research has succeeded in getting the optimal size and location of PV-DG as a solution to improve the performance of the distribution network on the Tipo system. The existence of a load growth factor can result in voltage drops on bus far from the power plant. By using NSGA-II, PV-DG of 8.2 MW spread over three buses can reduce power loss and improve the voltage profile. The PV-DG placement guarantees the voltage profile within tolerance limits even if there is an additional load until the third year.

The NSGA-II method shows better performance than the SGA method on the parameters of total output power, power loss reduction, and simulation processing time. The NSGA-II method shows that it requires faster convergence time processing than the SGA method with a greater reduction in system power loss.

REFERENCES

[1] Word Energy Outlook, 2014.
[2] Papermans, G., Driesen, j., Haeseldonckx, D., Belmans, R., & D, W., Distributed generation : definition, benefits and issue, Energy Policy, 33 (2005), 787-798.
[3] International Electrotechnical Commission, Efficient electrical energy transmission and distribution, Geneva, Switzerland: International Electrotechnical Commission, 2013.
[4] Fitriaty, P, & Zshen, Predicting energy generation from residential building attached Photovoltaic Cells in a tropical area using 3D modeling analysis, Journal of cleaner production, 195 (2017), 1422-1436.
[5] P.G.M Cormick, H.Suehrckeac, The effect of intermittent solar radiation on the performance of PV systems, Solar Energy, 171(2018), 667-674.
[6] M. C. V. Suresh, Edward J. Belwin, Optimal DG placement for benefit maximization in distribution networks by using Dragonfly algorithm, Renewables: Wind, Water, and Solar, 5 (2018), 1-8.
[7] Nirmalaj., V.Janamala and J. Rodrigues, Impact of Variable Distributed Generation on Distribution System Voltage Stability, International Conference on Data Science and Communication (IconDSC), 1-2 March 2019.
[8] J. A. Sa’ed, Q. Samara, S. Favuzza, G. Zizzo , Impact of integrating photovoltaic based DG on distribution network harmonics , IEEE International Conference on Environment and Electrical Engineering and 2017 IEEE Industrial and Commercial Power Systems Europe (EEEIC / I&CPS Europe), June (2017).
[9] Lasantha M., Tim L., Handbook of Distributed Generation, Distributed Solar-PV Generation: Impact on Voltage Control and Stability, (2017) ,317-342.
[10] Y. A. Rahman, S. Manjang, Yusran, and A. A. Ilham, Selection of Sensitive Buses using the Firefly Algorithm for Optimal Multiple Types of Distributed Generations Allocation, International Journal of Advanced Computer Science and Applications, 10 (2018), 316-322.
[11] Edmarcio A., Lina P., William M, Igor F., and Priscila R., Distributed Generation Allocation Using the Genetic Algorithm of Chu-Beasley and Sensitivity, Przegląd Elektrotechniczny, 8(2016).
[12] S.Daud, A.F.A.Kadir, C.K.Gan, A.Mohamed, Tamer Khatib, A comparison of heuristic optimization techniques for optimal placement and sizing of photovoltaic based distributed generation in a distribution system, Solar Energy, 140 (2016), 219-226.
[13] Rahman, Y.A., Manjang, S., Yusran, Ilham, A.A., An empirical metaheuristic assessment for solving of multi-type distributed generation allocation problem, Proceeding of IEEE International Seminar on Research of Information Technology and Intelligent Systems, ISRITI (2018).
[14] Arulraj, N. and Kumarappan.R., Optimal multiple installation of DG and capacitor for energy loss reduction and loadability enhancement in the radial distribution network using the hybrid WIPSO–GSA algorithm, International Journal of Ambient Energy, 41(2020), 129-141.
[15]Yiping, L., Gary G. Yen, Dunwei, G., A Multimodal Multiobjective Evolutionary Algorithm Using Two-Archive and Recombination Strategies, IEEE Transactions on Evolutionary Computation, 23 (2019 ).
[16] Partha P.Biswasa, P.N.Suganthana, Gehan A.J., Amaratunga, Decomposition based multi-objective evolutionary algorithm for windfarm layout optimization, Renewable Energy , 115 (2018), 326-337.
[17] B.Y.Quab J.J.Liangac, Y.S.Zhua, Z.Y.Wang, P.N.Suganthan, Economic emission dispatch problems with stochastic wind power using summation based multi-objective evolutionary algorithm, Information Sciences, 351(2016), 48-66.
[18] S. Kannan, S. Baskar, James D. McCalley, P. Murugan, Application of NSGA-II Algorithm to Generation Expansion Planning, IEEE Transactions on Power Systems, 24(2009 ).
[19] R. Y. Asmi, M. Salama, and I. A. Ahmad, Distributed generation’s integration planning involving growth load models by means of genetic algorithm, Archives of Electrical Engineering, 67(2018), No. 3, 667-682.
[20] Yusran, Yuli, A.R.,Indar C, Sri Mawar, S., Syafaruddin, Mesh grid power quality enhancement with synchronous distributed generation: optimal allocation planning using breeder genetic algorithm, Przegląd Elektrotechniczny, 1(2020), 82-86.

Authors: prof. dr. Ir. Salama Manjang, MT, Department of Electrical Engineering, University of Hasanuddin,Jalan Poros Malino km 6 Bontomarannu Gowa, Indonesia, E-mail: salamamanjang@unhas.ac.id; dr.Yuli Asmi Rahman,ST, M.Eng, Department of Electrical Engineering, University of Tadulako, Soekarno Hatta Km. 10, Tondo, Palu, Indonesia, E-mail: yuliasmi.rahman.81@gmail.com.

Source & Publisher Item Identifier: PRZEGLĄD ELEKTROTECHNICZNY, ISSN 0033-2097, R. 96 NR 9/2020. doi:10.15199/48.2020.09.12

Static Synchronous Compensator and Superconducting Fault Current Limiter for Power Transmission System Transient Stability Regulation Including Wind Generator

Published by Abdelkrim ZEBAR, Department of Electrical Engineering, Faculty of technology, Ferhat Abbas Setif 1 University, Setif, Algeria

Abstract. The increase of environmental pollution and the decreasing of the energy is extremely related to social progress. The resolution of the environment pollution depends on renewable energy such as wind energy system. Though, that system is required by transient and voltage stability issues with wind energy employing fixed-speed induction generators to be augmented with resistive type superconducting fault current limiter (SFCL) and dynamic compensation devices, such as a Static Synchronous Compensator units (STATCOM). The use of combined model based SFCL and STATCOM for promoting the transient and voltage stability of a multi-machine power system considering the fixed-speed induction generators is the primarily focus of this study. The proposed combined model functions top reserve simultaneously a flexible control of reactive power using STATCOM controller and to reduce fault current using superconducting technology based SFCL. The effectiveness of the proposed combined model is tested on the IEEE11-bus test system applied to the case of three-phase short circuit fault in one transmission line. A simulation results are presented in this document.

Streszczenie. Proponowany w artykule połączony model ma najwyższą rezerwę jednocześnie elastyczną kontrolę mocy biernej za pomocą kontrolera STATCOM i redukcję prądu zwarciowego za pomocą SFCL opartej na technologii nadprzewodnictwa. Skuteczność proponowanego modelu łączonego jest testowana na systemie testowym magistrali IEEE11 stosowanym w przypadku zwarcia trójfazowego zwarcia w jednej linii przesyłowej. Wyniki symulacji przedstawiono w tym artykule. (Statyczny synchroniczny kompensator i nadprzewodzący ogranicznik prądu dla systemu przesyłowego z generatorem wiatrowym)

Keywords: Distributed wind generation (DWG) , superconducting fault current limiter (SFCL), Static Synchronous Compensator (STATCOM), Transient stability.
Słowa kluczowe: Rozproszone system z z generatoirem wiatrowym (DWG), Nadprzewodnikowy ogranicznik prądu zwarciowego (SFCL), Statyczny synchroniczny kompensator (STATCOM), Stabilność przejściowa.

Introduction

Recently, the wind power generation included in the standard grid showed a further increase, in a significant way so, the electric utilities grid codes are forced to be revised, thus the reliability in systems with high wind energy diffusion will be guaranteed [1]. In power system stability studies the term transient stability usually refers to the ability of the synchronous machines to remain in synchronism during the brief period following large disturbances, such as severe lightning strikes, loss of heavily loaded transmission lines, loss of generation stations, or short circuits on buses [2]. To cope with the increasing demand for electric power, more and more FACTS devices are employed to improve the transmission capability of existing transmission facilities. As a result, the stability margin of power systems has decreased while the complexity of power systems has increased considerably. Thus, new techniques in power system control which can improve the dynamic performance and transient stability of power systems present an even more formidable challenge[3].

The STATCOM is one of the important FACTS devices and can be used for dynamic reactive power compensation of power systems to provide voltage support and stability improvement [4-5].

The use of Fault Current Limiters (SFCLs) is being evaluated as one element necessary to limit the fault current and enhance the power system transient stability [6]. A superconducting fault current limiter (SFCL) is a device with negligible impedance in normal operating conditions that reliably switches to a high impedance state in case of extra-current. Such a device is able to increase the short circuit power of an electric network and to contemporarily eliminate the hazard during the fault. It can be regarded as a key component for future electric power systems [7-8].

One of the requirements of transient stability analysis is to compute a transient stability index (TSI) for the contingencies, which is used to assess the stability of single contingency and furthermore rank the severity of different contingencies [9]. The Critical Clearance Time (CCT) of a fault is generally considered as the best measurement of severity of a contingency and thus widely used for ranking contingencies in accordance with their severity [10]. In this paper Critical Clearing Time (CCT) is employed as a transient stability index to evaluate test system. The Critical Clearing Time is defined as “the maximum time between the fault initiation and its clearing such that the power system is transiently stable”.

Many methods for transient stability analysis and assessment have been proposed and improved over the years, such as equal area criteria, numerical integration and Lyapunove method, in this study the numerical integration method is required in order to get the exact CCTs. The numerical integration method is the most reliable and accurate method for transient stability assessment [11].

Proposed work is aimed to investigate the potential influence of the combined application of superconducting fault current limiter (SFCL) and shunt FACTS Controller (STATCOM) for improving both transient stability and voltage regulation of the power system containing a distributed wind generation (DWG) based on conventional fixed speed induction generator. Moreover, the optimal location of the proposed coordinated controller (SFCL– STATCOM) is also analyzed. The effectiveness of the proposed combined model is tested on the IEEE 11-bus test system applied to the case of three-phase short circuit fault in one transmission line. Computer simulation results for system under study are presented and discussed.

Mathematical Model

This section gives a mathematical model for the power system network which includes modelling of synchronous generator, SFCLs and STATCOMs.

Synchronous Generator

With few characteristic assumptions, the synchronous generator can be modelled by the subsequent group of nonlinear equations [12-13]:

where: δ – rotor angle of the generator; ω – angular speed of the generator; ω₀ – rated generator rotor angle speed; H – inertia time constant; P_m – mechanical power; P_e = T_e ω₀ – – electrical power; D – damping coefficient; x_d – d-axis synchronous reactance; x_q – q-axis synchronous reactance ; x’_d – d-axis transient reactance ; x’_q – q-axis transient reactance ; E_d – d-axis transient voltage; E_q – q-axis transient voltage .

The electrical torque T_e as:

The equation of the stator are given by:

where: v_d, v_q : direct and quadrature axis stator terminal voltage components; V_t : generator terminal voltage; i_d, i_q : direct and quadrature axis stator current components.

Distributed wind generation

Distributed wind generation contain many wind turbines and their detailed modelling may be unaffordable due to computational burden. In order to reduce the dimensionality, aggregation techniques are used to obtain equivalent models. A proper equivalent model can be easily obtained for fixed-speed wind turbines where a one-to-one correspondence between wind speed and active power output exists. In this case, aggregation is performed by adding the mechanical power of each wind turbine and by using an equivalent squirrel cage induction generator (SCIG) which receives the total mechanical power [14-15].

A simplified transient model of a SCIG can be described by the following algebraic-differential equations [16]:

Here, x’ = (x_s + x_m x_r) / (x_m + x_r) is the transient reactance, R_s is the stator resistance which is assumed to be zero, x_r is the rotor reactance, xm is the magnetizing reactance, x = (x_s + x_m) is the rotor open circuit reactance, T₀‘ is transient open circuit time constant, Tm is the mechanical torque, S is the slip, T_e = (E_dr‘ I_ds – E_qr‘ I_qs) is the electrical torque, E_dr‘ and E_qr‘ are the direct and quadrature axis transient voltages respectively, I_ds and I_qs are the direct and quadrature axis currents respectively, and ω_s is the synchronous speed.

The DWG penetration level in the system is defined as [17]:

where P_DWG and P_CG are the amount of total active power generated by DWG and Centralized generation respectively.

Superconducting fault current limiter

Depending on the different superconducting materials and the operation principle the superconducting fault current limiters can be classified into different types [18]. In the resistive type the superconductor is directly connected in series to the line to be protected since in the inductive concept the superconductor is magnetically coupled into the line [19].

Fig.1. modified transmission line with SFCL

SFCL is a device that limits the fault current by generating impedance when a fault occurs. In addition, the limiting impedance generated to limit fault currents proves helpful in increasing generator output degraded by a fault, thus providing stabilization. as SFCLs installed in series with transmission lines can be just operated during the period from the fault occurrence to the fault clearing [7]. The equivalent circuit in π of the transmission line with SFCL is illustrated in Fig. 1, the associated equation for R_SFCL can be described by exponentially expressed as follows.

where, R_m is the expected maximum value of SFCL resistance in the normal state (R_m≈ 20 Ω), T_SC is the time constant of transition from the superconducting state to the normal state, which is assumed to be 1ms. The One phase of the resistance SFCL model is simulated in MATLAB/Simulink as shown in Fig. 2.

Fig.2. one phase of the resistance SFCL model

Static Synchronous Compensator

The Static Synchronous Compensator (STATCOM) is a shunt FACTS device that regulates the voltage of the ac bus to which it is connected. In the literature various STATCOM models have been developed and included within the load flow program, the optimal power flow and the transient stability analysis. A simplified STATCOM current injection model has been proposed in [20]. The STATCOM current i_sh is always kept in quadrature in relation to the bus voltage so that only reactive power is exchanged between the ac system and the STATCOM. The equivalent circuit and the control scheme are shown in Fig. 3 and undergoes the following differential equation:

where K_r is the regulator gain and T_r is the Regulator time constant.

The model is completed by the algebraic equation expressing the reactive power injected at the STATCOM node:

Fig.3. STATCOM circuit and control scheme

Simulation Results

To investigate the efficiency and the robustness of the proposed coordinated controller based SFCL and STATCOM on the power system transient stability enhancement in presence of distributed wind generation; the model is integrated in the IEEE benchmark four-machine two-area eleven bus test system in the case of three phase short circuit fault in the transmission line. DWG is connected to each of the load buses. The one line configuration is shown in Fig. 4.

Fig.4. on-line diagram of the Electrical Power System.

Technical data such as machine, voltage regulator, governor turbine, buses and branches information are taken from [21].

The transient stability is assessed by the criterion of relative rotor angles, using the time domain simulation method. The toolbox Sim Power Systems of MATLAB/SIMULINK software is used to carry out simulations studies.

Optimal Location Of SFCL–STATCOM

Optimal location and control of multi FACTS devices and multi SFCL is a vital and complex research area. In the literature many modern techniques and indices proposed for optimal location and control of multi FACTS devices. For secure operation of power systems, it is required to maintain an adequate voltage stability margin, not only under normal conditions, but, also, in contingency cases. In this study the voltage stability index using continuation power flow proposed for optimal location of STATCOM and SFCL.

From the continuation power flow results which are shown in the Fig. 5, the buses 5, 6, 7, 8, 9, 10 and 11 are the critical buses. Among these buses, bus 8 has the weakest voltage profile.

Fig.5. critical buses based on continuation power flow.

First, buses are classified based on three procedures: Procedure1: all buses are classified based on voltage stability index, the weak buses are identified based on voltage stability index, in this study, the bus 8 is considered as a candidate bus, the main role of the STATOM is to control voltage at this bus by exchanging reactive (capacitive or inductive) power with the network.

Procedure 2: Buses are classified based on the value of fault currents (three phase fault).

Procedure 3: Buses are classified based on the reactive power compensation witch consume by the DWG, the DWG will generate an active power, equal to the amount consumed by the load. However, in order to generate this necessary active power, the DWG need to consume reactive power from the network, the bus 9 is considered as the point of common coupling (PCC) where the WG is connected, the main role of the STATCOM is to compensate for this reactive power.

Impacts Of Combined SFCL–STATCOM Controller On Power System Transient Stability Enhancement

We have three logic Cases: Base case, which indicates the original system where there is no SFCL and STATCOM, in the system. Second case, with STATCOM at the weak bus (low voltage stability index) and SFCL at a bus which has high fault current. Third case, with STATCOM at the PCC and SFCL at another bus which has high fault current.

Case 1

A 3-phase fault is occurs at t = 1 second on line 7–8 near the bus 8 and it is cleared by opening the line at both ends. we consider a WG at bus 9 with a penetration level of 20 %. Generator 2 is the nearest generator to the fault location and therefore it has the most rotor speed deviation for this fault. The fault clearing time (FCT = 0.266 s) at first then (FCT = 0.300 s). Simulation results on the rotor angle differences of four generators without considering SFCL and STATCOM Controller are shown in Figs. 6.

Fig.6. relative rotor angles without SFCL–STATCOM.

It can be seen that the relative rotor angles are damped and consequently the system maintains its stability, but when the fault clearing time increased to 0.300 s, the relative angles (δ14, δ24 and δ34) increase indefinitely, so at this critical situation the system loses its stability.

Case 2

In order to maximize voltage stability index and to improve power system transient stability, STATCOM located at the weak bus (low voltage stability index) and the SFCL is placed in line 7–8 which has high fault current. The STATCOM will try to support the voltage by injecting reactive power on the line when the voltage is lower than the reference voltage. The first mentioned fault in the previous sub-section is applied again.

Time domain simulation performed at the cleared time 0.333 s, we can see from Fig. 8 the maximum relative rotor angles are (δ14 = 97°, δ24 = 91° and δ34 = 10°), the relative rotor angles (δ14, δ24 and δ34) are damped and therefore the system becomes stable compared to the first case (system unstable). The current on line 7–8 with SFCL is shown in Fig. 7.

In Fig. 9 it can be also seen that the system response with the SFCL-STATCOM is better than that with the with only STATCOM (Fig. 8) in the sense of the settling time is reduced. The critical clearing time is enhanced to a new value (0.355 s).

Fig.8. Relative rotor angles with only STATCOM.

Fig.9. Relative rotor angles considering SFCL– STATCOM.

Case 3

In case 2 the SFCL is placed in the line 7–8 which has high fault current and the STATCOM located at the weak bus however in this case the STATCOM is placed at the PCC. The multipurpose becomes: reduce the current in line 7–8 (high fault current) and maximize dynamically voltage stability index, in this case the STATCOM compensate the reactive power consumed by the DWG and the fault current reduced by SFCL enhance the performance of the STATCOM (reduce saturation problem) dynamically during fault, and alternatively the required size of STATCOM will be reduced (economic aspect).

Fig.10. Relative rotor angles considering SFCL– STATCOM

As a result, the reactive powers delivered by generating units reduce. Compared to the two others cases, the critical clearing time is improved. The SFCL is placed in line 7–8. The first mentioned fault in the sub-section (case 1) is applied again. The fault is cleared after 0.427 s.

In Fig. 10, It can be seen that the maximum relative rotor angles are (δ14 = 20°, δ24 = 18° and δ34 = 3°), the relative rotor angles (δ14, δ24 and δ34) are damped and therefore the system becomes stable compared to the first and second cases (system unstable). It can be also seen that the system response with the STATCOM at the PCC is better than that with the STATCOM at the weak bus in the sense of the settling time is reduced. The critical clearing time is enhanced to a new value (0.482 s). It is important to conclude that integration of shunt FACTS compensator (STATCOM) in coordination with SFCL in suitable location may help the system to improve the transient stability. Table 1 shows the values of margin stability (CCT) obtained corresponding to different cases.

Table 1. Margin stability (CCT) for different cases

Conclusion

Modern electric power systems’ operation, control and stability have been heavily affected by The rising penetration of energy sources renewal, increasing demands, restricted resources, and deregulated electricity markets power systems. In this study, the multi-machine power system transient stability improvement contains a large DWG via superconducting fault current limiter (SFCL) and shunt FACTS Controller (STATCOM) when applied through coordinated application was discussed and investigated. Simulation results performed on the IEEE benchmarked four-machine two-area test system in presence of distributed wind generation considering three phase short circuit clearly indicate that proposed combined controller placed at suitable locations can be used as an effective mean capable to enhance the margin stability and extend the critical clearing time in a multi-machine power system.

Future research will focus on investigation the effect of combined application of superconducting fault current limiter and other transient stability improvement FACTS devices in the presence of the distributed generation by considering the optimal value of SFCL and location of this hybrid Controller.

REFERENCES

[1] M. F. Farias, M. G. Cendoya, P. E. Battaiotto, Wind Farms in Weak Grids Enhancement of Ride-Through Capability Using Custom Power Systems, IEEE/PES Transmission and Distribution Conference and Exposition Latin America, pp. 1-5, 2008.
[2] A. Karami and S.Z. Esmaili, “Transient stability assessment of power systems described with detailed models using neural networks”, International Journal of Electrical Power and Energy Systems 45 (2013) 279–292
[3] L. Cong, Y. Wang, D.J. Hill.Transient stability and voltage regulation enhancement via coordinated control of generator excitation and SVC. Journal of Electrical Power and Energy Systems 27 (2005) 121–130
[4] R. You, M. H. Nehrir & D. A. Pierre, ” Controller Design for SVC and TCSC to Enhance Damping of Power System Oscillations”, Electric Power Components and Systems, Volume 35, Issue 8, pages 871-884, May 2007
[5] A.E. Hammad, “Analysis of power system stability enhancement by static var compensators”, IEEE Trans. PS 1 (4) (1986) 222–227.
[6] AL. Bettiol, A. Souza, JL. Todesco, JR. Tesch., “Estimation of critical clearing times using neural network. In: IEEE bologna powertech conference”, Bologna, Italy; 2003.
[7] M. Sjostrom, R.Cherkaoui and B.Dutoit, “Enhancement of power system transient stability using superconducting fault current limiters”, IEEE Trans. Applied Superconductivity, vol.9, no.2, pp.1328-1330, 1999.
[8] M. MAJKA, J. KOZAK, The Coreless Superconducting Fault Current Limiter 15 kV 140 A. Przegląd Elektrotechniczny, 92 (2016), no. 7, 38- 41
[9] K. Phorang, M. Leelajindakraireak, and Y. Mizutani, “Damping improvement of oscillation in power system by fuzzy logic based SVC stabilizer” Asia Pacific. IEEE/PES Transmission and Distribution Conference and Exhibition 2002, Vol.3, Oct. 2002, pp.1542–1547
[10] R. Ebrahimpour, E. K. Abharian, S. Moussavi, Z. & A. A Motie Birjandi, “Transient stability assessment of a power system by mixture of experts” International Journal of Engineering (IJE) Volume (4): Issue (1). pp. 93–104
[11] P. Kundur, “Definition and classification of power system stability IEEE/CIGRE joint task force on stability terms and definitions” IEEE Transactions on Power Systems, vol. 19, no.2, p. 1387-1401, May 2004
[12] N. ABU-TABAK, “Stabilité dynamique des systèmes électriques multi-machines: modélisation, commande, observation et simulation”, Doctoral Thesis, University of Lyon, 2008.
[13] I. GRICHE, S. MESSALTI , K. SAOUDI, Parallel Fuzzy Logic and PI Controller for Transient Stability and Voltage Regulation of Power System Including Wind Turbine. Przegląd Elektrotechniczny, 95 (2019), no. 9, 51- 56
[14] H. Arnaldo, Wind farm model for power system stability analysis, Doctoral Thesis, University of Illinois at Urbana- Champaign, 2010.
[15] A.Zebar, A. Hamouda and K. Zehar, “Impact of the location of fuzzy controlled static var compensator on the power system transient stability improvement in presence of distributed wind generation”, Rev. Roum. Sci. Techn. – Électrotechn. et Énerg., 60, 4, p. 426–436, Bucarest, 2015.
[16] Roy, Naruttam Kumar, M. J. Hossain, and H. R. Pota, Voltage profile improvement for distributed wind generation using DSTATCOM, IEEE Power and Energy Society General Meeting, 2011.
[17] M. Reza, PH. Schavemaker, JG. Slootweg, WL. Kling, L. Van Der Sluis, Impacts of distributed generation penetration levels on power systems transient stability, IEEE Power Engineering Society General Meeting, 2004.
[18] M. Noe, M. Steurer, “High-temperature superconductor fault current limiters: concepts, applications, and development status”, SUST 20 (2007).
[19] S. Nemdili and S. Belkhiat, “Electrothermal modeling of coated conductor for a resistive superconducting fault-current limiter”, J. Supercond. Nov. Magn. (2012). doi:10.1007/s10948-012- 1895-4
[20] F. Milano, Power System Modelling and Scripting, London: Springer, Aug. 2010.
[21] M. Klein, G. Rogers, Moorty and P. Kundur: “Analytical investigation of factors influencing PSS performance” , IEEE Trans. on Energy Conversion, vol. 7, no 3, pp. 382-390, September 1992.

Author: Dr. Abdelkrim. ZEBAR is with the Department of Electrical Engineering, Faculty of technology, and Ferhat Abbas Setif1 University, Setif, Algeria (corresponding author to provide phone: +213-776-87-72-46; e-mail: karimzebar@univ-setif.dz).

Source & Publisher Item Identifier: PRZEGLĄD ELEKTROTECHNICZNY, ISSN 0033-2097, R. 96 NR 11/2020. doi:10.15199/48.2020.11.33

Wind Power Forecasting Based on Meteorological Data Using Neural Networks

Published by 1. Yuriy SAYENKO¹, 2. Ryszard PAWEŁEK², 3. Vadym LIUBARTSEV¹,
Pryazovskyi State Technical University (1), Lodz University of Technology, Institute of Electrical Power Engineering (2) ORCID. 1. 0000-0001-9729-4700, 2. 0000-0002-1023-8210, 3. 0000-0003-1243-9101

Abstract. The growing share of renewable energy sources in the structure of energy systems causes many problems related to the correct operation of the grid. This impact is most evident in low-voltage grids to which many low-power prosumer solar and wind installations are connected. For the correct management and, consequently, the economic operation of power systems, the most accurate forecast of electricity consumption and generation in grids with different voltage levels is needed. Conventional generation devices have stable production values and can be regulated within wide limits, while the production of electricity from renewable sources, by wind farms in particular, depends on external weather conditions and requires a more careful approach to its forecasting. The aim of the article is to present a method of forecasting the power generated by wind turbines based on publicly available meteorological data. The presented forecasting method uses the theory of neural networks.

Streszczenie. Rosnący udział odnawialnych źródeł energii w strukturze systemów energetycznych, powoduje wiele problemów związanych z poprawną pracą sieci. Oddziaływanie to jest najbardziej widoczne w sieciach niskiego napięcia, do których przyłączonych jest wiele fotowoltaicznych i wiatrowych instalacji prosumenckich małej mocy. Dla poprawnego zarządzania i w konsekwencji ekonomicznej pracy systemów elektroenergetycznych potrzebna jest możliwie dokładna prognoza zużycia i wytwarzania energii elektrycznej w sieciach o różnych poziomach napięcia. Konwencjonalne urządzenia wytwórcze mają stabilne wartości wytwarzania i mogą być regulowane w szerokich granicach, natomiast produkcja energii elektrycznej ze źródeł odnawialnych, a w szczególności przez elektrownie wiatrowe, zależy od zewnętrznych czynników atmosferycznych i wymaga staranniejszego podejścia do jej prognozowania. Celem artykułu jest przedstawienie metody prognozowania mocy generowanej przez turbiny wiatrowe w oparciu o publicznie dostępne dane meteorologiczne. W prezentowanej metodzie prognozowania wykorzystano teorię sieci neuronowych. (Prognozowanie energii wytwarzanej przez źródła wiatrowe na podstawie danych meteorologicznych z wykorzystaniem sieci neuronowych).

Keywords: renewable sources, wind farms, forecasting, neural networks, modelling
Słowa kluczowe: źródła odnawialne, farmy wiatrowe, prognozowanie, sieci neuronowe, modelowanie

Introduction

Nowadays, renewable energy sources occupy a large part of the generation structure. The advantages of using renewable types of energy sources are as follows: reduced CO₂ emissions and the reduced consumption of hydrocarbon fuels, but they also have a significant drawback, i.e. complex management and dispatching of such power systems due to the instability of electricity production [1 – 4].

At the same time, if large wind-powered generating plants have the equipment and information resources to predict generation depending on weather conditions, small wind-powered generating plants (the prosumers) do not have such capabilities and the usual meteorological weather forecast is the only available data source. The power system operator must have the data forecast on possible short-term and long-term generation from the wind-powered generating plants in order to maintain the most economical operating mode of the power system with distributed generation, without reducing its reliability indicators.

The power production of the wind-turbine depends on many variables, primarily such as wind speed, wind direction, temperature, air pressure, etc., which must be taken into account in the forecasting model.

Currently, the greatest attention is devoted to the issues of forecasting electricity production from large wind-powered generating plants that have the meteorological equipment and the necessary information resources to transfer the electricity production plan to the grid operator [5 – 9].

But at the same time, there is an issue of forecasting electricity production from the prosumers of small wind turbines as part of an electric power system with distributed generation, because they are also capable of having a strong effect on the power flow and the operating mode of the power system.

The aim of the article is to develop a method of predicting the power generation from wind turbines that is based on public meteorological data.

Research method

Nowadays, neural networks are one of the most progressive and accurate methods of forecasting and assessing various processes in power engineering. The feasibility of their use is confirmed by previous studies [10, 11], in which their main advantages are determined: high accuracy, versatility and flexibility.

A neuron is an information–processing unit that is fundamental to the operation of a neural network. The block diagram of Figure 1 shows the model of a neuron which forms the basis for designing a large family of neural networks studied in the next chapters. Here, we identify three basic elements of the neural model [12]:

– A set of synapses, or connecting links, each of which is characterized by a weight or strength of its own. Specifically, signal x_j at the input of synapse j connected to neuron k is multiplied by the synaptic weight w_kj. It is important to take note of the manner in which the subscripts of the synaptic weight w_kj are written. The first subscript in w_kjrefers to the neuron in question and the second subscript refers to the input end of the synapse to which the weight refers. Unlike the weight of a synapse in the brain, the synaptic weight of an artificial neuron may lie in a range that includes negative as well as positive values.

– An adder for summing the input signals, weighted by the respective synaptic strengths of the neuron; the operations described here constitute a linear combiner.

– An activation function for limiting the amplitude of the output of a neuron. The activation function is also referred to as a squashing function, in that it squashes (limits) the permissible amplitude range of the output signal to some finite value.

The neural model of Figure 1 also includes an externally applied bias, denoted by b_k. Bias b_k has the effect of increasing or lowering the net input of the activation function, depending on whether it is positive or negative.

In mathematical terms, we may describe neuron k depicted in Figure 1 by writing the pair of equations [11]:

and

where: x₁, x₂, …, x_m – input signals; w_k1, w_k2, …, w_km – respective synaptic weights of neuron k; v_k – linear combiner output due to the input signals; b_k – bias; φ(·) – activation function; y_k – output signal of the neuron.

Source data

Data on electricity production taken from the website of the Australian energy market was used as a source of preliminary data for the training of the neural network. The data is a time series for the generation of electricity by various wind-powered generating plants in Australia with a discreteness of 5 minutes for 2018 (about 105 000 values). The data obtained from 195 MW wind farm in Portland was used in the study. This choice is due to the rather close (16 km) location of the weather station at the Portland airport (Fig. 2), and it was data freely available for research. In preparation, the data was reduced to the discreteness of 1 hour.

Fig.2. The relative position of the wind farm and the meteorological station in Portland (Australia)

The archive of the meteorological data obtained from the meteorological station located at the Portland airport [14] includes data on time, temperature, air pressure, wind speed and direction, as well as pressure trends (Fig. 3).

The time data is presented in fractions of the whole value of the day (24 hours is 1) for the ease of use in the neural network. The wind direction data for the training of neural networks is presented as 0.01 part of the wind rose angle (for example, north wind – 0, east – 0.9, south – 1.8). The rest of the data was accepted unchanged.

Fig.3. The block diagram of the source data for the neural network training

Creation and training of neuron networks

To predict the power generation via meteorological data, a neural network with the following variables was established (the selection of the variables was based on the fact that they provide high speed and quality prediction selected for different configurations):

– Number of layers – 2;
– Function of the first layer activation: Log–sigmoid;
– Function of the second layer activation: Linear;
– Training algorithm – Bayesian regularization back propagation.

The number of neurons N was selected empirically [9], based on the following ratio:

where N_IN– number of inputs (in the analysed case N_IN = 6).

A simplified model of a neuron network created in the MATLAB program using the Neural Network Toolbox package [15] is presented in Figure 4.

Fig.4. A simplified structure of a neuron network for forecasting electricity generation by wind farms

The performance graph (Fig.5) indicates the iteration at which the validation performance reached a minimum. The training continued for 6 more iterations before it was stopped.

In the analysed case, the figure does not indicate any major problems with the training. The validation and test curves are very similar. If the test curve had increased significantly before the validation curve increased, it would have been possible for some overfitting to have occurred.

The next step in validating the network is to create a regression plot which shows the relationship between the outputs of the network and the targets (Figure 6). If the training were perfect, the network outputs and the targets would be exactly equal, but the relationship is rarely perfect in practice.

The dashed line in each plot represents the perfect result – outputs = targets. The solid line represents the best fit linear regression line between outputs and targets. The R value is an indication of the relationship between the outputs and targets. If R = 1, this indicates that there is an exact linear relationship between outputs and targets. If R is close to zero, there is no linear relationship between outputs and targets. For this example, the training data indicates a good fit. The validation and test results also show R values close to 0.9. The scatter plot is helpful in showing that certain data points have poor fits.

Fig.6. Graphs of training, validation and testing of the neuron network

Figure 7 shows the graph with an example of real (blue) and predicted values obtained as a result of the neural network operation. Forecasting was carried out during training on different volumes of data sets (1500–5500 values).

The forecasting was carried out 100 hours ahead, then the values were compared with the actual 100 values that were previously received from the Australian grid operator.

Fig.7. An example of a forecast of electricity generation in comparison with real values

Table 1 shows the average generation power forecast errors for wind-powered generating plants. The following equation was used to calculate them:

where: Err is the error value, MW; P_pred is the predicted value of the generation power, MW; P_real is the real value of the generation power, MW.

Next, for each element in the error vector, the error value is calculated as a percentage:

where: ErrPer – is the error value, %; P – is the installed capacity of the wind-powered generating plant, MW.

To get the final result of the average vector error as a percentage, the following equation was used:

where MeanErrPer – is the average error value, %; N – is the number of predicted values of electricity generation.

The comparison of the average errors when using different ranges of data sets is given in Table 1:

Table 1. Comparison of the values of the mean forecast errors

Conclusion

The use of neural networks to predict the generation of electricity from wind-powered generating plants, including private ones of low power, provides an accurate forecast (1- 2% average error).

It is possible to obtain an accurate forecast using conventional meteorological weather forecast data just for private small wind farms without meteorological stations (data on temperature, air pressure, and most importantly, wind speed and direction) as the initial data.

When the forecasting uses neural networks, it is necessary to take into account the requirements for the quality of the initial data and its volume. For example, when using a sample of data of 1500 values for training a neural network, an average error of 14.3% was obtained, and when a sample of 2500 values was used, it was 1.15% (Table 1).

At the same time, the error value increased to 4.91% with too large a sample (5500 values), which is caused by over-generalization for the time series.

Therefore, the objective is to find a compromise for the sample size for training the network. This forecast can also be used by operators of distributed generation power systems to assess the impact of weather conditions on the power generation capacity by prosumers and more accurately balance power flows to achieve the most economical power system operation mode.

REFERENCES

[1] Leithon J., Werner S., Koivunen V., Cost-aware renewable energy management: Centralized vs. distributed generation, Renewable Energy, pp. 1164-1179, 2019 doi:10.1016/j.renene.2019.09.077
[2] Wang C., Wu J., Ekanayak J., Smart Electricity Distribution Networks, Boca Raton: Taylor & Francis Group, 2017.
[3] Kariniotakis G., Waldl I.H-P., Marti I., Giebel G., Nielsen T.S., Tambke J., Usaola J., Dierich F., Bocquet A., Virlot S., Next generation forecasting tools for the optimal management of wind generation, 2006 International Conference on Probabilistic Methods Applied to Power Systems, Stockholm, 2006, pp. 1-6, doi: 10.1109/PMAPS.2006.360238.
[4] Shahnia F., Arefi A., Ledwich G., Electric Distribution Network Planning, Power Systems, 2018, doi:10.1007/978- 981-10-7056-3
[5] Zhang Y., Dong J, Least Squares-based Optimal Reconciliation Method for Hierarchical Forecasts of Wind Power Generation”, IEEE Transactions on Power Systems, pp. 1–10, 2018, doi:10.1109/tpwrs.2018.2868175
[6] Yatiyana E., Rajakaruna S., Ghosh A., Wind speed and direction forecasting for wind power generation using ARIMA model, 2017 Australasian Universities Power Engineering Conference (AUPEC), Melbourne, VIC, 2017, pp. 1-6, doi: 10.1109/AUPEC.2017.8282494.
[7] Tambke J., von Bremen L., Barthelmie R., Palomares A.M., Ranchin T., Juban J., Kariniotakis G.N., Brownsword R.A., Waldl H.P., Short-term Forecasting of Offshore Wind Farm Production – Developments of the Anemos Project”, 2006 European Wind Energy Conference, EWEC, Feb 2006, Athènes, Greece, 13 p.
[8] Wan C., Xu Z., Pinson P., Dong Z.Y., Wong K.P., Probabilistic Forecasting of Wind Power Generation Using Extreme Learning Machine, IEEE Transactions on Power Systems, vol. 29, no. 3, pp. 1033-1044, May 2014, doi: 10.1109/TPWRS.2013.2287871.
[9] Foley A.M., Leahy P.G., Marvuglia A., McKeogh, E.J., Current methods and advances in forecasting of wind power generation”, Renewable Energy, vol. 37(1), pp. 1–8, 2012, doi: 10.1016/j.renene.2011.05.033
[10] Sayenko Y., Sychenko V., Liubartsev V., Development of Methods for Optimizing Reactive Power Modes Based on Neural Network Technologies, 2019 IEEE 6th International Conference on Energy Smart Systems (ESS), Kyiv, Ukraine, 2019, pp. 98-103, doi: 10.1109/ESS.2019.8764220.
[11] Sayenko Y., Baranenko T., Liubartsev V., Forecasting of Electricity Generation by Solar Panels Using Neural Networks with Incomplete Initial Data, 2020 IEEE 4th International Conference on Intelligent Energy and Power Systems (IEPS), Istanbul, Turkey, 2020, pp. 140-143, doi: 10.1109/IEPS 51250.2020.9263085.
[12] Haykin S.O., Neural Networks and Learning Machines, 3rd Edition. Ontario, Canada: McMaster University, 2009.
[13] AEMO Energy Generation Data, Australian Energy Market Operator. [Online]. Available: https://anero.id/energy/data
[14] Weather archive in Portland (airport), Australia, [Online].available: https://rp5.ua/Weather_archive_in_Portland_(airport)_Australia
[15] Hudson Beale M., Hagan M.T., Demuth H.B., Neural Network Toolbox™. Users’s guide, The MathWorks, Inc., 2017.

Authors: prof., Ph.D., D.Sc., Eng., Yuriy Sayenko, Pryazovskyi State Technical University, Department of Industrial Electrical Power Supply, Ukraine, 87555, Mariupol, 7 Universytets’ka, E-mail: sayenko_y_l@pstu.edu; Dr Ryszard Pawełek, Lodz University of Technology, Institute of Electrical Power Engineering, 18/22 Stefanowskiego str., 90- 924 Lodz, E-mail: ryszard.pawelek@p.lodz.pl; Ph.D. student., Vadym Liubartsev, Pryazovskyi State Technical University, Department of Industrial Electrical Power Supply, Ukraine, 87555, Mariupol, 7 Universytets’ka, E-mail: lubartsevvadim@gmail.com.

Source & Publisher Item Identifier: PRZEGLĄD ELEKTROTECHNICZNY, ISSN 0033-2097, R. 97 NR 11/2021. doi:10.15199/48.2021.11.39

An Application of Facility Location Problem for Electricity Supply Cost Minimization at the Stage of Preliminary High Voltage Network Development Planning

Published by Piotr KAPLER, Warsaw University of Technology, Faculty of Electrical Engineering, Power Engineering Institute

Abstract. The article deals with minimizing the costs of electricity supply at the preliminary network planning stage. It presents selected information on costs in power engineering. The facility location problem and Mixed-Integer Linear Programming has been also described. A computational example of the application of this method is presented to solve the problem of minimizing the costs of electricity supply between high voltage / 110 kV power stations and 110 kV urban stations

Streszczenie. Artykuł dotyczy minimalizacji kosztów dostaw energii elektrycznej na etapie wstępnego planowania sieci. Przedstawiono w nim wybrane informacje dotyczące kosztów w elektroenergetyce. Opisano zagadnienie lokalizacji obiektów oraz metodę optymalizacji mieszanej całkowitoliczbowej liniowej. Zaprezentowano przykład obliczeniowy zastosowania tych metod do rozwiązania problemu zminimalizowania kosztów dostaw energii elektrycznej pomiędzy stacjami NN / 110 kV a stacjami miejskimi 110 kV. (Zastosowanie problemu lokalizacji obiektów do minimalizacji kosztów dostaw energii elektrycznej na etapie wstępnego planowania rozwoju sieci wysokiego napięcia).

Keywords: costs of electricity, costs minimizing, network development planning, mixed integer linear programming
Słowa kluczowe: koszty energii elektrycznej, minimalizacja kosztów, planowanie rozwoju sieci, optymalizacja mieszana.

Introduction

Electric Power System (EPS) consists of many interconnected elements whose purpose is the generation, transmission, distribution and usage of electricity. Simultaneously, it is a very important and strategic element of any area. Therefore, it is necessary that the decisions concerning power system taken at each stage of its design, operation or development, are optimal and reasonable. Thanks to this approach it should be possible to rationally manage energy sources, good operation of components as well as sustainable use of electricity by consumers both now and in the future.

Due to its complexity and large scale electric power system is a challenge for many optimization processes. Optimization tasks in this case can be both linear and nonlinear and have a large computational dimension. They are also accompanied by a number of different constraints, that must be maintained in order to model the correct operation of the power system. Often, the variables involved in the optimization process can be either integer or non-integer. Also some problems are NP-hard type. Optimization is proposed in each of the power system sectors: from generation to end consumers of electricity. Moreover, optimization may apply to both planned and already existing power infrastructure. Usually time horizon of optimization can be middle (2-5 years) or long term (10-20 years). The network planning process itself sometime takes into account the conditions of risk [1], uncertainty [2] or a probabilistic approach [3]. Factors encouraging the development are: an increase in power demand, the location of a new generation sources, access to fuels or increasing reliability of supplies [4].

This article is devoted to the minimization of electricity supply costs by application of facility location problem. Described situation occurs at the stage of preliminary high voltage (HV) network planning with usage of mixed integer linear optimization. The initial part of article presents selected issues about costs in power engineering. Also a description of the optimization method used is provided. The calculation part of article presents an example of minimizing the costs of electricity supply between potentially planned high voltage / 110 kV power stations and 110 kV urban stations. The final part of the article contains a summary and conclusions.

Selected issues concerning the costs of electricity supply

The operation of the power system is based on continuous coverage of the variable power demand of end users. Electricity, by its nature, cannot be stored easily and cheaply at present. At the same time, its production and transmission are associated with significant costs. These premises encourage rational management of electricity.

The supply of electricity is possible by building an appropriate infrastructure consisting of lines and power stations. The financial outlays for these investments are very large. In addition, it is also necessary to earn money on the provision of transmission services and collect funds for the ongoing operation and renovation works. In the future, investments in further network development may be required. Money can be obtained due to the difference between the total revenues that arose from the sale of goods or services and the total costs that had to be incurred in producing those goods and services.

In the case of power lines, expenditure on their construction is related with the length of the planned line, land charges and the use of specific technical solutions (for example types of poles or wires). It is estimated that the financial expenditure on the construction of a cable line is a multiple of 4 to 14 times the expenditure on the construction of an overhead line with the same rated voltage and length [5]. The operating costs result from the need to conduct qualified service work. Modifying the connections later on is cheaper and easier for overhead lines than for cable lines.

For power stations, investment outlays are related to the area where station is to be located, with the equipment used (for example circuit breakers or switch disconnectors) and the type of transformers. Operating costs will depend on the station layout used.

In addition, the transmission of electricity is accompanied by energy losses that must be paid for. Their values will depend on the length of the connections. They can be divided into load losses and idle losses.

Various methods are used in the electricity sector to determine fees for the provision of transmission services. The most popular include: the postage stamp method, the contract path method, the MW-mile method or method of tracing power flows. Each of the above methods give different outcomes as a result of the adopted assumptions and simplifications in their operation. The mentioned methods can be divided according to the way costs are treated. There are embedded cost methods and marginal cost methods. The first subgroup includes the postage stamp, contract paths and MW-mile methods. The second one is power flow tracing. Embedded costs deal with the power system as a whole and usage costs are borne by all users [6]. The marginal cost methods determine the value of the unit price to cover variable costs.

Linear mixed integer optimization and facility location problem

Linear programming has been used to solve optimization problems for many years. It is a special case of mathematical programming. Increased interest in this subject occurred at the turn of the 1950s and 1960s. Despite the fact that it is not new issue, nowadays it sill has a great application due to the development of computers and the application of possibilities for modern technical problems.

The idea behind this tool is to build a mathematical model that best reflects the characteristics of a real object as much as possible. Such a model is a record in the form of appropriate mathematical equations. The solution of the equations is the answer to the given decision problem. This approach helps making the final decision because building a mathematical model is cheap, simple and safe and does not require manipulating real objects.

In general, the linear programming problem can be defined in matrix form as (1):

subject to conditions (2) and (3):

where: z – objective function, c’ – matrix of objective function coefficients, A – matrix of limiting conditions coefficients, x – matrix of variables (unknowns), b – matrix of limiting values of the right side of the inequality (2). Formula (3) is a non-negative condition.

Each model consists of three basic elements. These are: the objective function (also known as the criterion function), decision variables and set of limiting conditions. The objective function is describing the goal to be pursued in the optimization task. This function is being minimized or maximized. It has a linear form. Decision variables are a description of resources related to the modelled issue. They are non-negative. The constraint conditions reflect the limits in the modelled process. They are described as inequalities with left sides containing linear functions.

Despite the simplicity of implementation, linear programming in its classic form is not always suitable for solving most optimization problems. It may turn out that in the description of the problem it will be necessary to include binary values (0 or 1) representing for example the existence or non-existence of a given object depending on certain conditions. In order to solve such problems, linear mixed integer optimization (MILP – Mixed Integer Linear Programming) is used. This approach applies if some (but not all) of the variables are binary.

In general, MILP problems can be defined identically using formula (1) and (2). However, condition (3) is replaced with condition (4).

where some variables are integer type and the rest are binary (0 or 1).

This type of approach is particularly useful for solving optimization problems in the power engineering industry, for example: the switching on or off generation units [7], designing distribution networks [8] or multistage planning of the development of the transmission network [9].

Optimization problems solved with the use of MILP must meet the following simplifying assumption: the share of each of the variables in the objective function is proportional to this variable and, at the same time, independent of the values of other variables in the task. Additionally, the values of constraints and decision variables must be known before performing optimization calculations. Sometimes, however, these values are not known exactly because they were derived from a rough estimate. Their influence on the optimization result can be obtained by applying a sensitivity analysis.

The facility location problem is well-known in optimization theory. In general, it deals with selecting the best locations of facilities that will cover the given customers demand while minimizing the costs of transportation. The input data are: set of candidate facility locations F = {1,…,m}, set of customers C = {1,…,n}, a cost function and a distance function. Moreover, the cost function may contain not only information about the cost of delivery but also the cost of opening a given facility. The result of calculation is the assignment of all customers to all opened facilities. Optimization models with facility locational problem can also be extended with additional constraints like: availability of fuels, raw materials or cost of land. Further description of this problem can be found in [10,11].

There are many variations of this task. In presented solution the mentioned problem has a MILP form – opening (or not) a new facility is a binary value. The objective is to minimize the total costs which can be divided to costs of opening new high voltage / 110 kV power station and costs of supply (which are proportional to the distance between power station and customers).

Formulation of the optimization problem

The optimization problem concerned the issue of minimizing the cost of electricity supply at the stage of preliminary high voltage network planning. The solution was obtained with an application of facility locational problem for power system planning. In the presented case facilities were 4 potentially planned high voltage / 110 kV power stations while customers were 9 urban 110 kV power stations. Every customer station can be theoretically supplied by each of planned HV / 110 kV station. The aim of the task was to find such connections between power stations so that the cost of electricity supply was minimal while simultaneously covering the required demand for active power.

All 110 kV power stations were located at different distances from supplying stations so the delivery cost was also dissimilar. It was assumed that the transmission of electricity to all customers can take place without violating any technical limitations like voltage drops or long-term current carrying capacity. Furthermore, it was also assumed that the total active power in supply stations was equal to the demand of consumers (Case 1) or supply exceeds that demand (Case 2).

The objective function (minimizing the cost of electricity supply) was defined as follows (5):

subject to conditions (6), (7) and (8):

where: C – objective function constituting the cost of electricity supply subject to minimization, S – total number of power stations i, b_s – the price of the building of the i-th power station, k_su – binary value related to the use of the station (0 – station not used, 1 – station used), e_s – cost of i-th station exploatation, c – delivery cost matrix for each station-customer pair, i – the amount of electricity transmitted from the i-th station to j-th customer, k_l – power loss cost, C_s – total number of customers j, a_s – active power available at the station, d_s – demand for active power at given customer.

In the objective function (5) it was taken into account that the cost of electricity supply from the HV / 110 kV to 110 kV stations will depend on the distance, active power demand, fixed cost related to the operation of a given high voltage station and the value of the charge for transmission power losses. Table 1 contains the values of the cost of delivery parameters for all potential station-customer pairs. They were created by applying the MW-mile method in which the delivery cost is defined by the formula (9) [12]:

where: C_dⁱ – cost Energy supply, P_jⁱ – active power flow in the j-th line for the i-th transaction, D_j– length of the j-th line, R_j– required unit renewal value per line length.

The unit value of the renovation was 1.5% of the power line construction cost. The average price of building 1 km of high voltage lines was determined on the basis of data from [13]. The cost of construction of the high voltage / 110 kV power station was estimated at PLN 150 million (around 37 million euros). The costs of maintenance and repairs of each power station were assumed for 4% of capital expenditure [14].

Table 1. Energy supply costs form the MW-mile method, all values in PLN.

At the preliminary planning stage, it was assumed that all possible station-to-customer connections were taken into account. Some of these connections may later be rejected in further stages of network planning if it turns out that they are less optimal than others.

Solving the optimization problem

The task of minimizing the cost of electricity supply was solved in two cases. Case 1 assumed that the sum of active power in all power stations was equal to the sum of active power demand for all customers. In case 2, it was assumed that the planned total active power in all stations is greater than the customers demand. Table 3 presents the active powers in the stations for case 1 and case 2.

The solution of the problem was presented in the form station-customer pairs, for which, as a result of optimization calculations, active power values were allocated to cover the total demand at the customer (110 kV station) with the lowest possible cost of delivering and meeting additional conditions – expressed by the equations (6)-(8).

Table 4 shows binary values (1 or 0) corresponding to the necessity to build (or not) a given high voltage / 110 kV station in each of the cases. Table 5 contains the final summary of the results for both analysed cases. Where the connection of a given station with customer would not be optimal, the value of 0 appears. It is worth noting that connections with the value of 0 can be both in case 1 and case 2. This means that regardless of the values of available active power the given connection is not optimal.

Table 2. Demanded active powers, in MW

Table 3. Active powers in high voltage power stations, in MW

Based on the presented optimization results, the following values of the objective function (electricity supply cost) were obtained: for case 1 – PLN 650 244 922, and for case 2 – PLN 493 955 302. The result for case 2 is smaller because the calculations show that building of high voltage station 3 is redundant. Consequently, each Station 3 – customer pair has a 0 value in Table 5. Hence the result for case 2 is taken as an optimal for the considered planning problem.

Table 4. Binary values corresponding to the need to build a given power station in each case.

Table 5. Results of the optimization task solution, active powers in MW.

Summary and conclusions

The article presents an optimization task consisting in minimizing the costs of electricity supply at the stage of preliminary high voltage network planning. According to the conducted research, case 2 turned out to be more favourable than case 1. The objective function for case 2 is around 75% of the cost for case 1. However, the difference in cost values is PLN 156 289 620. This is more than the assumed construction price of the one high voltage / 110 kV power station.

The presented example relates to the situation when it is planned to transmit power form high voltage / 110 kV power stations to 110 kV stations located in a large urban agglomeration. In this case, the knowledge from the above analysis may be useful in further planning stages of the power network development for example, to the construction of load flow and short-circuit network models. The mentioned method may also be effective for making decisions concerning medium and low voltage network development. Due to the use of binary values, it can also be used to decide on the closing of existing power stations.

The results of the performed analysis depend on the input data used. It is necessary to choose the method by which the charges for transmitting electricity will be calculated. It is assumed that the price criteria used in the planning process should also take into account the problem of renewal of elements [15]. The second important factor is knowing the prices for the construction of the power line. The value of this price will strongly depend on many factors like type of line (overhead or cable) or the way it is located.

The advantage of the presented methodology using the mixed-integer approach is its efficiency, scalability and simplicity. The model used is universal and after introducing minor modifications it can be adopted to other optimization problems in the power engineering industry. Even if the obtained solution does not turn out to be fully satisfactory, it may constitute a starting point for the search of a new result. The disadvantage of an inaccurately conducted analysis may be both underinvestment and overinvestment in the transmission network [16]. Also, using the mixed-integer approach may be difficult due to possibly large number of binary (0,1) decision values.

When solving the problems of planning the power network development it is necessary to conduct multivariant analyses of various possible technical solutions. On this basis, it will be possible to further conclude whether it is worth implementing the intended investment or not. In the process of building decision models it is convenient to use, for example, linear programming. In addition, knowledge about the necessary investment outlays is also essential.

REFERENCES

[1] Marzecki J., Planowanie rozwoju miejskich Rozdzielczych Punktów Zasilających (RPZ) w warunkach ryzyka, Przegląd Elektrotechniczny, 90 (2014), nr. 2, 234-237
[2] Marzecki J., Planowanie rozwoju miejskich stacji 110 kV/SN w warunkach niepewności, Przegląd Elektrotechniczny, 96 (2020), nr. 1, 23-26
[3] Kubek P., Przygrodzki M., Wybrane aspekty wykorzystania elementów probabilistycznych w planowaniu rozwoju sieci przesyłowej, Przegląd Elektrotechniczny, 94 (2018), nr. 12, 108-111
[4] Gonen T., Electrical Power Transmission System Engineering: Analysis and Design, CRC Press, 2014
[5] Underground vs. Overhead: Power Line Installation-Cost Comparison and Mitigation, https://www.powergrid.com/2013/02/01/underground-vs-overhead-power – line-installation-cost-comparison/#gref (accessed 24.04.2020)
[6] Saganek D., Koszty wykorzystania elementów SEE w powiązaniu ze stanami pracy SEE – właściwości podejścia opartego na idei wpływu, Przegląd Elektrotechniczny, 91 (2015), nr. 5, 159-165
[7] Kamiński J., Kaszyński P., Mirowski T., Szurlej A., Krótkoterminowy matematyczny model systemu wytwarzania energii elektrycznej dla warunków Polski, Rynek Energii, czerwiec 2014
[8] Turkay B., Distribution System Planning Using Mixed Integer Programming, Elektrik, 6 (1998), nr. 1, 37-48
[9] Zhang H., Vittal V., Heydt G. T., J. Quintero, A mixed-integer linear programming approach for multi-stage securityconstrained transmission expansion planning, IEEE Transactions on Power Systems, 27 (2012), nr. 2, 1125-1133
[10] Dantzig G. B., Linear Programming and Extensions, Princeton University Press, Princeton, New Jersey, 1991
[11] E. Castillo, A. J. Conejo, P. Pedregal, R. Garcia, N. Alguacil, Building and Solving Mathematical Programming Models in Engineering and Science, John Willey and Sons, 2002
[12] Song Y-H. (edit), Modern Optimisation Techniques in Power Systems, Springer Science + Business Media, 1999
[13] Ciupak S., Linie WN na słupach kratowych vs linie WN na słupach rurowych, Konferencja Naukowo-Techniczna Elektroenergetyczne linie napowietrzne i kablowe wysokich i najwyższych napięć, Wisła 18-19.10.2017
[14] E. Dyka, I. Mróz-Radłowska, Ekonomia w Energetyce. Wybrane zagadnienia, Politechnika Łódzka, Łódź 2014.
[15] Einhorn M., Siddiqi R. (edit), Electricity Transmission Pricing and Technology, Kluwer Academic Publishers, 1996
[16] Dołęga W. Aspekty rynkowe planowania rozwoju sieciowej infrastruktury elektroenergetycznej, Energetyka, 8 (2015)

Author: Piotr Kapler, Ph. D., Warsaw University of Technology, Faculty of Electrical Engineering, Power Engineering Institute, Koszykowa 75, 00-662 Warsaw, Poland. E-mail: piotr.kapler@ien.pw.edu.pl

Source & Publisher Item Identifier: PRZEGLĄD ELEKTROTECHNICZNY, ISSN 0033-2097, R. 96 NR 11/2020. doi:10.15199/48.2020.11.35

Methods of Calculating Solar Insolation for the Assessment of Energy Efficiency of Solar Power Plants

Published by Sergey CHIZHMA, Artyom ZAKHAROV, Immanuel Kant Baltic Federal University (Kaliningrad, Russia)

Abstract. The authors explore a number of methods of calculating insolation on a tilted arbitrarily oriented surface and choose the most efficient one for conducting the performance and efficiency analysis of solar power plants. The validity of the method is tested using the statistical data of the Kaliningrad region. The authors describe an algorithm for calculating insolation in the Matlab software programme, based on the following data: the position of the Sun, the average monthly climatic characteristics of the region, beam, scattered and reflected insolation as well as the proportion of these types of insolation on a tilted surface. The authors validate the results of their calculations by comparing them with statistical averages.

Streszczenie. Analizowano metody obliczania nasłonecznienia dowolnie pochylonej powierzchni. Walidację metody przeprowadzono na podstawie danych statystycznych regionu Kaliningrad. Obliczenia przeprowadzono na podstawie znajomości pozycji słońca, średniej miesięćżnej informacji klimatycznej, średniego nasłonecznienia. (Metoda obliczania nasłonecznienia dowolnie pochylonej powierzchni)

Keywords: alternative energy, insolation, energy, efficiency
Słowa kluczowe: energia nasłonecznienia, ogniwa fotowoltaiczne

Introduction

The construction of solar power plants requires assessment of energy potential of a chosen territory. In most cases, design solutions for solar power plants are based on long-term average measurements that include the average daily solar energy (kWh/(m²·day)), average daily wind speed (m/s), etc. The main sources of data used for estimating solar energy potential are insolation maps, NASA data [1] and the data provided by the World Radiation Data Centre [2] for actinometric stations of the World Meteorological Organization (WMO) network, etc. [3]. There is an extensive body of literature on hybrid energy solutions, types of renewable energy, their classification and composition. Reviews on different types of solar energy solutions are presented in [4-6].

There is a need for a more accurate assessment of the efficiency of solar power plants. This task requires more objective data reflecting the impact of a number of factors associated with the location and the position of solar panels, annual variations in the length of daylight hours, and cloudiness. In our research, we propose to move from the average energy parameters, which are, in essence, quantitative characteristics of insolation, to instantaneous values of insolation incident on a tilted surface in order to determine daily fluctuations of insolation. We also propose a method for assessing the efficiency and performance of solar power plants based on this approach.

The aim of the research

The main aim of our research is to test a method of calculating insolation for obtaining more accurate data on its instantaneous values based on geographical location and average climatic parameters. We hold that the obtained data can be used for modeling and assessing the energy generation performance of solar power plants [7] and the optimization of their control algorithms. We tested the proposed method using input data for the Kaliningrad region and validated the results of the calculation.

To achieve the aim, we analysed a number of calculation methods, which are based on the algorithm for determining the position of the Sun for different geographic locations and seasons (Section 1). Various computational methods for estimating insolation on the Earth’s surface are described in the works of Liu and Jordan [8-10], KolaresPereira and Rabl [11], and Gueymard [12-13].

We have chosen one calculation method that meets the objectives of our research (See Section 2) and used it for the determination of solar insolation on a tilted surface (Section 3). The selected model has been applied to obtain insolation data for the Kaliningrad region for one year. Power plants based on renewable energy sources can be simulated with obtained data, as shown in [14-15]. The results of our calculations have been verified (Section 4).

1. Methods of calculating the position of the Sun

The existing methods of estimating instantaneous values of insolation require the analysis of the position of the Sun above the Earth’s surface [16]. The annual change in solar activity is associated with two main factors: the annual change in the solar declination angle δ_s, and the annual change in the distance between the Earth and the Sun. The axial tilt of the Earth’s orbit is 23.45°. Thus, the declination angle δ_s (solar declination) varies in the range from -23.45° to 23.45° during the year and can be calculated using the Cooper formula [17]:

where n is the ordinal number of the day in the year.

To simplify the calculation of the position of the Sun above the Earth’s surface, the Earth is considered as a stationary object. Figure 1 shows that the observation point is assumed as the origin of the coordinates. In this system, the Sun is a moving object and rotates around the centre of the Earth with an angular velocity of 15° per hour. The zenith point (or local noon) is taken as the starting point for the revolution of the Sun. The value characterizing the movement of the Sun around the Earth is called the solar hour angle h_S. It is calculated using the formula:

The two values, the angular height of the Sun α and the azimuth angle of the Sun α_s are used to calculate the spherical coordinates of the position of the Sun above the Earth’s surface:

where L is the geographical latitude.

Fig.1. Determination of the solar hour angle h_s, the angle of declination of the Sun δ_s, the angular height of the Sun α and the azimuth angle α_s

Solar radiation on the Earth’s surface depends on the changes in the distance between the Earth and the Sun during the year: on the day of the winter solstice, December 21, it is 1,471·10¹¹ m and on the day of the summer solstice, June 21, it constitutes 1,521·10¹¹ m. This variability is a function of the elliptical shape of the Earth’s orbit, which is taken into account when translating the so-called solar time (the zenith time is 12.00) into local time using the following formula:

where ET is the equation of time, taking into account the variability of the Earth’s rotation speed around the Sun, l_st is the standard time meridian, l_local is the local longitude. ET is determined using the formula given below:

2. The analysis of the methods of calculating insolation on the Earth’s surface

The calculation of insolation on the surface of the Earth is a difficult task since it requires taking into account a variety of factors. There are many empirical or semiempirical approaches to the problem within which a number of simplifying assumptions are made. Each of the assumptions is based on the calculation of three main types of insolation on the Earth’s surface: beam insolation (I_b.c), scattered or diffused insolation (I_d.c) and reflected insolation (I_r.c):

To account for the absorbing effect of clouds, a monthly clearness index (K_T) is introduced. It can be calculated using the formula:

where H_h [kWh/(m²·day)] is the monthly average daily radiation on a horizontal surface, and H_o.h [kWh/(m²·day)] is the monthly average daily extraterrestrial radiation [16].

This index was first proposed in the works of Liu and Jordan [8–10]. The method of calculating radiation, considered in the framework of this approach, is applied to monthly average values. In our research, we propose to calculate the instantaneous values of insolation during the day. There are a number of methods that allow researchers to move from average values of solar radiation to instantaneous ones:

• the Kolares-Pereira and Rabla method developed in 1979 (CPR-method) [11],

• the Kolares-Pereira and Rabla method, modified by Gyeumard in 1986 (CPRG-method) [12],

• the Daily Integration (DI) Model, developed by Gyeumard in 2000 [13].

Based on the data presented by Christian A. Gueymard in [13], the daily integration (DI) model is more accurate than other methods listed above. Therefore, in this paper, the calculation of insolation is done using this method. The daily integration model is based on the calculation of r_d and r_t, which change during daylight hours. They demonstrate the dependence of instantaneous diffused insolation and instantaneous insolation on a horizontal surface on the total daily scattered insolation and total daily insolation on a horizontal surface, respectively.

The coefficient r_d is defined as follows:

where h_ssis the solar light angle corresponding to the sunset, hereinafter expressed in radians. The method of calculating h_ss is described in [16].

According to the daily integration model, the r_tcoefficient is calculated as follows:

where the absorption of radiation in the atmosphere is taken into account using the ratio a₂/a₁:

The coefficients A(h_ss) and B(h_ss) are defined as follows:

Missing values for the calculation are defined as indicated below:

Thus, the instantaneous value of insolation on the horizontal surface I_h [Wh/m²] and the instantaneous value of diffused insolation I_d [Wh/m²] can be calculated in the following way:

Then, beam insolation on a horizontal surface, I_d.h [Wh/m²] is calculated according to the formula:

3. Calculation of solar radiation on a tilted surface

The values of insolation obtained using the daily integration method are to be converted into insolation values for solar panels oriented and tilted in a particular way. The position of solar panels is set at the following angles (Fig. 2): the tilt of the solar panels β and the azimuth angle of the solar panels α_w. The incident angle of sunlight on the surface of the solar panels i is determined as follows:

Using simple transformations, we determine the fraction of beam insolation on the tilted surface I_b.c:

Fig.2. Calculation of insolation on a tilted surface

Scattered insolation (I_d.c) on the surface of solar panels is estimated as follows:

Reflected insolation (I_r.c) is calculated using the formula:

where ρ = 0.2 for the surface not covered with snow, and 0.8 for the snow-covered surface.

Thus, the insolation on the tilted surface is determined as follows:

4. Calculation of insolation for the Kaliningrad region

The calculation method described in the sections above was used in the MATLAB software programme for obtaining the annual insolation data for the Kaliningrad region (54°43′N, 20°30′E). The algorithm (Fig. 3) is based on NASA resource data [1], where H_h is average monthly total radiation on a horizontal surface per day, H_d is average monthly total scattered radiation, and K_T is monthly clearness index:

Fig.3. Block diagram of the algorithm for calculating insolation

Other input data for the algorithm are the latitude, the longitude of the chosen terrain, the start and finish dates of the formation of the calculated insolation values, as well as parameters characterising the location of the solar panels. The algorithm performs the calculation of the ordinal numbers of the days in the year (taking into account leap years) for which the calculation of solar insolation is done. NASA data are monthly averaged, so before starting the calculation, the change in monthly average indicators is interpolated, and their values for each day are determined. The variable, linearly dependent on time, is the hour angle of the Sun h_s, for which every minute values are calculated. After all the data have been formed, the MATLAB software programme calculates the total insolation on a tilted surface, thus forming the data array I_c.

Figure 4 shows the calculated insolation on the surface of a solar panel having a 30° tilt angle on December 21, June 21, and March 1.

The reliability of the calculation method and the data obtained have been assessed. To perform the assessment we analysed the calculated total monthly radiation incident on a horizontal as well as on a tilted surface. The results are shown in Fig. 5, where H_h.calcuis the calculated average monthly total radiation, H_h is the average monthly total radiation according to NASA [1], and H_c is the monthly average total radiation.

Fig.5. Estimation of the data reliability

The graph shows that the margin of error of the proposed method is minimal. It allows us to use the obtained data for the performance and efficiency analysis of solar power plants.

Conclusions

In our research, a number of methods for calculating instantaneous values of insolation have been analysed. The methods are based on the calculation of the coordinates of the position of the Sun (the angular height and the azimuth angle), as well as the average climatic parameters of the season.

Our analysis of the calculation methods has shown that the daily integration model (DI) is the most effective one. On the basis of this model, we proposed an algorithm for calculating insolation using the MATLAB software programme. This allowed us to determine instantaneous values of insolation for a wide variety of temporal and geographic conditions. The algorithm can be used to form an array of objective data on the instantaneous values of insolation on the tilted surface of solar panels, as a combination of beam, diffused and reflected insolation.

We validated the obtained data and identified the margin of error of the proposed method, which turned out to be insignificant. The proposed method of calculation can be used to simulate the operation of power plants working on renewable energy sources, evaluate their efficiency and to conduct further research on optimizing control algorithms for autonomous solar power plants.

REFERENCES

[1] https://power.larc.nasa.gov/data-access-viewer/.
[2] World Radiation Data Centre. Available at: http://wrdc.mgo.rssi.ru/ (accessed 14 March 2017)
[3] Popel’ O.S., Frid S.E., Kiseleva S.V., Kolomiec Ju.G., Lisickaja N.V. Klimaticheskie dannye dlja vozobnovljaemoj jenergetiki Rossii (Baza klimaticheskih dannyh). M.: Izd-vo MFTI. 2012
[4] Chauhan A., Saini R.P. A review on Integrated Renewable Energy System based power generation for stand alone applications: Configurations, storage options, sizing methodologies and control // Renewable and Sustainable Energy Reviews. – 2014. –V. 38. – Р. 99–120.
[5] Shivarama K.K., Sathish K.K. A review on hybrid renewable energy systems // Renewable and Sustainable Energy Reviews. – 2015. – V. 52. – Р. 907–916.
[6] Badwawi R.A., Abusara M., Mallick T. A Review of Hybrid Solar PV and Wind Energy System // Smart Science. – 2015. – V. 3 (3). – Р. 127–138.
[7] Obuhov S.G., Plotnikov I.A. Imitacionnaja model’ rezhimov raboty avtonomnoj fotojelektricheskoj stancii s uchetom real’nyh uslovij jekspluatacii. Izvestija Tomskogo politehnicheskogo universiteta. Inzhiniring georesursov. 2017. T. 328. № 6. 38–51
[8] Liu, B.Y.H., R.C. Jordan. 1960. The interrelationship and characteristic distribution of direct, diffuse and total solar radiation. Sol. Energy 4: 1–19.
[9] Liu, B.Y.H., R.C. Jordan. 1961a. Daily insolation on surfaces titled toward the equator. Trans. ASHRAE 67: 526–541.
[10] Liu, B.Y.H., R.C. Jordan. 1961b. Daily insolation on surface tilted toward the equator. Trans. ASHRAE 3(10): 53–59
[11] Collares-Pereira M, Rabl A. The average distribution of solar radiation: Correlations between diffuse and hemispherical and between daily and hourly insolation values. Solar Energy 1979; 22: 155-164.
[12] Gueymard C. Monthly averages of the daily effective optical air mass and solar related angles for horizontal or inclined surfaces. J Solar Energy Eng Trans ASME, 1986
[13] Gueymard C. “Prediction and Performance Assessment of Mean Hourly Global Radiation” Solar Energy, Vol. 68, No. 3, 2000, pp. 285-303. doi:10.1016/S0038-092X(99)00070-5
[14] Mirosław Mazur , Janusz Partyka , Tomasz Marcewicz. Analysis of the use of a hybrid power system of renewable wind and photovoltaic energy in residential buildings. PRZEGLĄD ELEKTROTECHNICZNY, ISSN 0033-2097, R. 92 NR 8/2016. 113-116.
[15] Kazimierz Buczek, Wiesława Malska, Sebastian Penar. Use of PSIM software for modelling a small solar power station. PRZEGLĄD ELEKTROTECHNICZNY, 08/2011 Page no. 42.
[16] D. Yogi Goswami. Principles of solar engineering. Third Edition. CRC Press. Taylor & Francis Group 2015.
[17] Duffie J.A., Beckman W.A. Solar Engineering of Thermal Processes. Hoboken, New Jersey, John Wiley & Sons, Inc., 2013.

Authors: Prof. Sergey N. Chizhma, the Institute of Physics, Mathematics and Information Technologies, Immanuel Kant Baltic Federal University, Kaliningrad Russia. E-mail: chisn@yandex.ru, Artyom Zakharov, PhD student, the Institute of Physics, Mathematics and Information Technologies, Immanuel Kant Baltic Federal University, Kaliningrad, Russia. E-mail: AIZakharov@kantiana.ru.

Source & Publisher Item Identifier: PRZEGLĄD ELEKTROTECHNICZNY, ISSN 0033-2097, R. 96 NR 11/2020. doi:10.15199/48.2020.11.26

The Problem of Determining the Coefficient of Flicker in Accordance to Normative Regulations

Published by Marta BĄTKIEWICZ-PANTUŁA, Wrocław University of Science and Technology, Institute of Electrical Power Engineering

Abstract. The article presents an analysis of the selected parameter of the power quality. The analysis was done in accordance with the applicable normative regulations. The results presented in the article were made on the basis of actual measurements of the electricity quality parameters. Two different objects were analyzed: a group of industrial receivers and a network cooperating with renewable energy sources.

Streszczenie. W artykule zaprezentowano analizę wybranego parametru jakości energii elektrycznej. Analiza została przeprowadzona zgodnie z obowiązującymi przepisami normatywnymi. Zamieszczone w artykule wyniki zostały wykonane na podstawie rzeczywistych pomiarów parametrów jakości energii elektrycznej. Analizie zostały poddane dwa różne obiekty jakim była grupa odbiorników przemysłowych oraz sieć współpracująca z odnawialnymi źródłami energii. (Problematyka wyznaczania współczynnika migotania światła zgodnie z aktualnymi przepisami normatywnymi).

Keywords: flicker, regulation of the minister, updating of the EN 50160 standard, industrial customer of electricity, renewable power plants
Słowa kluczowe: współczynnik migotania światła, rozporządzenie ministra, aktualizacja normy EN 50160, odbiorcy przemysłowi, energia odnawialna

Introduction

The phenomenon of voltage fluctuations appeared with the beginning of the existence of electricity distribution systems. More attention to this phenomenon have been given with the increase in the number of receivers and their installed capacity. In order to better understand the disorders and the effects they cause, in many countries they started to fight this problem through numerous researches aiming at a thorough understanding of the phenomenon itself and its measuring assessment.

Until recently, voltage fluctuations in the power supply or the load terminals, were described by specifying the maximum size of the change in value of the effective voltage.

The energy of voltage fluctuations and the spectral density of power of voltage fluctuations were also used, (also called the energy spectrum of voltage fluctuations). The duration of the fluctuations and the interval between changes in voltage was also included in the assessments.

Parameters according to which the quality of the supply voltage should be assessed are included in the Regulation of the Minister of Economy [1], which is a valid legal act and in relation to recipients supplied from low, medium and high voltage public power networks in PN-EN 50160 [2,3]. An important group of documents in this area is the multipart PN-EN 61000-xx standard with the general title “Electromagnetic compatibility” [4,5], which defines the acceptable levels of environmental and operational interference that are required to ensure the correct operation of electrical devices connected to the network and defines the measurement methods and methods for determining these disturbances.

The phenomenon of the flicker is a voltage changes that have the character of regular fluctuations, which occurs for a long time, which can cause changes in luminous flux generated by electric light sources. This phenomenon has an adverse effect on the comfort and concentration of people working in such conditions, and its measure is the flickering index, the permissible levels are specified in the Regulation [1] and the Standard [2]. The indicator of flicker is determined for all voltage levels in the power network [2], both for the highest voltages (connection groups I and II), as well as medium and low voltage (connection groups III – V).

The determination of acceptable levels of the flicker indicator also in the case of medium and high voltages, i.e. those at which the light source is not directly supplied, indicates the high importance of this parameter and the possibility of propagating interference at various voltage levels across the entire distribution network. Connecting the devices of significant powers to the network, which may cause the phenomenon of flickering of light (welding machines, often switched on and off motors, arc furnaces) should be preceded by an analysis aimed at checking whether the flickering effect caused by this device will fall within the permissible range [2]. In the analysis of this important parameters are: the nominal power of the device, the short-circuit power of the network and the nature of voltage changes caused by the operation of the device.

Measurement of voltage fluctuations is carried out in order to assess the compliance of existing levels of the phenomenon with the relevant standards, as well as to determine the emission level of a given receiver and compare it with the limit values in the standards. There are two basic methods of measurement:

• based on a quantitative evaluation of the phenomenon based on a temporary change in the effective value or voltage envelope,

• based on indirect measurement – measurement of the phenomenon of flickering of light which is a direct result of voltage fluctuations.

The phenomenon of flickering of light, known in the world literature under the name “flicker”, is one of the parameters for assessing the power quality. It depends on the instability of the perception of human vision, but this instability is caused by a light stimulus whose luminance or spectral distribution is subject to changes in time due to fluctuations in the voltage supplying the light source.

The rules for calculating the flicker ratio are specified in the standard [6]. This factor consists of two elements:

• short-term flicker factor P_st, determined for the observation time of 10 minutes, according to the dependence:

in which P_0.1 percentiles; P₁; P₃; P₁₀ and P₅₀ are flicker levels exceeded by 0.1; 1; 3; 10 and 50% of the observation time. The index s in the above dependence indicates that should used the smoothed values.

• long-term flicker factor P_lt determined for the observation time of 2 hours and calculated using the next 12 P_st factors for this observation time, according to the formula:

The primary document in the process of assessing the power quality is the Regulation of the Minister of Economy of 4 May 2007 on detailed conditions for the operation of the power system with the last update in 2008. The subordinate document is the PN-EN 50160 standard: Voltage characteristics of electricity supplied by public electricity networks which until now has been updated three times, i.e. in 2008, 2010 and 2015.

The Regulation of the Minister of Economy and the PNEN 50160 standard to be updated in 2010 were in line with the power quality requirements for the assessment.

Table I presents the comparison of the requirements for the parameter of the power quality which is the flicker. The list has been presented for nn and SN networks.

Table 1. Permissible limits of the flicker factor

The first column indicates the requirements set by the Regulation of the Minister of Economy in the second and third column requirements to be met according to PN-EN 50160 in 2010 and from 2015.

Analysis Of Measurements

The basis for assessing the power supply conditions is the Regulation of the Minister of Economy of 4 May 2007. on detailed conditions of the power system operation. The article presents the assessment of the parameters of power quality in accordance with the Regulation of the Minister [1] and the PN-EN 50160 [2] standard and its update [3]. According to the normative assumptions, a representative period was selected which was a normalized time basis for the presented runs and determined power quality indicators. In a representative measurement period if disturbances were noted, i.e. dips, power interruptions, according to [3], no such measurement results should be considered. The article uses only those measurement results that were considered significant from the point of view of this work.

Figure 1 shows the general schematic diagram of the measurement system. The main element of the measurement system was the Fluke 1760 recorder.

Fig.1. General schematic diagram of the measuring system

The actual measurements presented in the article have been divided into two groups. In the first group of industrial receivers supplied from the power grid was presented two examples of circuits, which are, the welding circuit, and installed arc furnace circuit. In the second group are the working with renewable sources medium voltage energy networks, also present in two cases, with the network hydroelectric plant and a wind power plant.

The first of the discussed cases are groups of industrial receivers supplied from power grid.

Figure 2 can be seen the values of the long-term and short-term flicker indicators on the L1 phase, obtained during measurements for welding circuit. The full observation time of the recorded parameters was one week. In the representative measurement period, no disturbances were detected that would affect the analyzed process. The indicators that are calculated in accordance with the methods described in the standard [4] are shown on below waveforms.

Fig.2. Flicker on example of L1 phase: Short-term indicator P_st (thin); P_lt long term indicator (bold)

It can be observed that for the analysed case, the power quality parameter does not exceed the permissible limits contained in the documents [1-3]. The value of the flicker indicator (long-term) is a maximum of 0.85 which is in accordance with the required ministerial regulation [1], standard [2] which refers to 95% of observation time and norm [3] for 100% of observation time. Permissible limits for the short-term coefficient were not specified in the Minister’s Regulation [1] and the Standard [2]. Coefficient P_st = 1.4 meets the requirements of the updated standard [3] despite exceeding the required value 1.2. The value was exceeded in less than 5% of cases, which corresponds to requirement 1.2 for 95% of the measurement data set.

The next analyzed case is presented in Figure 3. The recorded values of long-term and short-term flicker on the L1 phase example, were presented, they were collected during measurements for the circuit where arc furnace was installed. The complete observation time of the recorded parameters was four days. In the representative measurement period, no disturbances were detected that would affect the analyzed process. The presented waveforms refer to flicker indicators was calculated in accordance with the methods described in the standard [4].

Fig.3. Flicker on example of L1 phase: Short-term indicator P_st (thin); P_lt long term indicator (bold)

It can be observed that for the analysed case, the power quality parameter does not exceed the permissible limits contained in the documents [1-3]. The value of the long-term flicker factor is a maximum of 0.13 which is in accordance with the required ministerial regulation [1], standard [2] which refers to 95% of observation time and norm [3] for 100% of observation time. Permissible limits for the short-term flicker factor were not specified in the Minister’s Regulation [1] and the Standard [2]. The P_st = 0.16 coefficient meets the requirements of the updated standard [3].

The second of the discussed cases are medium voltage network working with renewable energy sources The first case where renewable sources are working of with medium voltage network is presented in Figure 4. Presented are the collected values of long-term and short-term flickers on the example of phase L1, obtained during measurements for a hydroelectric power plant. The full observation time of the recorded parameters was one week. In the representative measurement period, no disturbances were detected that would affect the analyzed process. The presented waveforms refer to the flicker factors calculated in accord with the methods from the standard [4].

It can be observed that for the analyzed case the power quality parameter exceeds the acceptable limits contained in the documents [1-3]. Value of the long-term flicker is a maximum of 1.8 which is in accordance with the required ministerial regulation [1], standard [2] which refers to 95% of observation time. Standard [3] requires fulfillment of the parameter for 100% observation time, so exceeding the permissible value of 1 to 1.8 results in failure to meet the acceptable limit. The coefficient of short-term flicker in accordance with the regulation of the minister [1] and the norm [2] does not have admissible limits. The P_st = 4.0 coefficient meets the requirements of the updated standard [3] despite exceeding the required value 1.2. The value was exceeded in less than 5% of cases, which corresponds to requirement 1.2 for 95% of the measurement data set. In the analyzed case, a solution should be considered to improve the power quality parameters.

Fig.4. Flicker on example of L1 phase: Short-term indicator P_st (thin); P_lt long term indicator (bold)

The next case is presented in Figure 5. Presented are the collected values of long-term and short-term flicker factors on the example of phase L1, obtained during measurements of a wind farm. The full observation time of the recorded parameters was one week. In the representative measurement period, no disturbances were detected that would affect the analyzed process. The presented waveforms refer to the flicker factors determined in accord with the methods from the standard [4].

Fig.5. Flicker on example of L1 phase: Short-term indicator P_st (thin); P_lt long term indicator (bold)

It can be observed that for the analysed case, the power quality parameter does not exceed the permissible limits contained in the documents [1-3]. The value of the long-term flicker indicator is a maximum of 0.8 which is in accordance with the required ministerial regulation [1], standard [2] which refers to 95% of observation time and norm [3] for 100% of observation time. Permissible limits for the short-term coefficient were not specified in the Minister’s Regulation [1] and the Standard [2]. The P_st = 1.8 coefficient meets the requirements of the updated standard [3] despite exceeding the required value 1.2. The value was exceeded in less than 5% of cases, which corresponds to requirement 1.2 for 95% of the measurement data set.

Summary

The tightening of the normative provisions that followed the introduction of the PN-EN 50160 [3] standard was aimed at improving the parameters of the power quality.

The change in the power quality parameters concerned the setting of a limit for the short-term flicker factor, which until the update was undefined. At the same time, the requirements set for the long-term factor flickers were tightened. Additional information that introduces the update of the standard is defining the normal working time. The concept of normal working time refers to the state of the systems without any disturbance. If during the analyzed measurement period there are events, they should be excluded from further analysis.

Analyzing the examples presented in the article, it can be concluded that the tightening of the limits permitted for the flicker factors did not significantly affect the power quality for individual or industrial consumers. The tightening of limits has affected the network cooperating with renewable sources. Renewable power plants are characterized by variability of work which is associated with a large irregularity of parameters. The tightening of allowable limits has had a particularly strong effect on rapidly variable and irregular burdens. It can be concluded that for customers where device start-ups will occur frequently or there will be large parameter changes, the introduction of a limit for the short-term flicker factor will result in non-compliance with the normative requirements.

The problem with tightening the permissible limits for power quality parameters occurs when according to the superior document, i.e. the minister’s regulation [1], the power quality parameters meet the requirements and subordinate document, i.e. the PN-EN 50160 [3] standard, the power quality parameters already conditions do not meet. Until the update in 2015, such a problem did not exist because the subordinate document defined the same limits as in the parent document.

REFERENCES

[1] Regulation of the Minister of Economy dated. May 4, 2007 on detailed conditions of functioning of the power system (Journal of Laws of 29 May 2007, item 623).
[2] EN 50160:2010. Voltage characteristics of electricity supplied by public electricity networks.
[3] EN 50160:2015. Voltage characteristics of electricity supplied by public electricity networks.
[4] PN-EN 61000-4-15: 2011 Electromagnetic Compatibility (EMC) – Test and Measurement Methods – Flicker Meter – Functional and design specifications.
[5] PN-EN 61000-4-30:2015 Electromagnetic compatibility (EMC)-Part 4-30: Testing and measurement techniques – Power quality measurement methods
[6] PN-EN 61000-3-3:2013, Electromagnetic compatibility – Permissible levels – Limiting voltage fluctuations and flicker caused by receivers with rated current < or = 16 A in low voltage supply networks.
[7] Guide to Quality of Electric Supply for Industrial Installations, Part 5, Flicker and Voltage Fluctuations, “Power Quality” Working Group WG2, 2000
[8] A.Klajn, M.Bątkiwicz-Pantuła, Application Note – Standard EN50160: Voltage characteristics of electricity supplied by public electricity network,. European Copper Institute, 2013.

Authors: dr inż. Marta Bątkiewicz-Pantuła, Wrocław University of Science and Technology, Faculty of Electrical Engineering, Institute of Electrical Power Engineering, 27 Wybrzeże Wyspiańskiego St, 50-370 Wrocław, Poland. E-mail: marta.batkiewiczpantula@pwr.edu.pl

Source & Publisher Item Identifier: PRZEGLĄD ELEKTROTECHNICZNY, ISSN 0033-2097, R. 96 NR 1/2020. doi:10.15199/48.2020.01.13

Analysis of Transient Waveforms in a Power System at Asymmetrical Short-Circuits

Published by Piotr PRUSKI, Stefan PASZEK, Silesian University of Technology

Abstract. In the paper, there is presented the analysis of disturbance waveforms of a synchronous generator operating in a single-machine power system consisting of a generating unit connected by a high-voltage power line to a bus. There are considered disturbances in the form of symmetrical and asymmetrical short-circuits in different places of the transmission line. The current-voltage equations of the power line and the bus are written for phase components, which allows for easy modeling of various asymmetries.

Streszczenie. W artykule analizowano przebiegi zakłóceniowe generatora synchronicznego pracującego w jednomaszynowym systemie elektroenergetycznym składającym się z zespołu wytwórczego, połączonego linią energetyczną wysokiego napięcia z siecią sztywną. Uwzględniono zakłócenia w postaci zwarć symetrycznych i niesymetrycznych w różnych miejscach linii przesyłowej. Równania prądowo-napięciowe linii energetycznej i sieci sztywnej zapisano dla składowych fazowych, co pozwala na łatwe modelowanie różnych asymetrii. (Analiza przebiegów nieustalonych w systemie elektroenergetycznym przy zwarciach niesymetrycznych).

Keywords: power system; generating unit; asymmetrical short-circuits; transient states.
Słowa kluczowe: system elektroenergetyczny, zespół wytwórczy, zwarcia niesymetryczne, stany nieustalone.

Introduction

Short-circuits are the most common and severe kind of faults occurring in the power system (PS). Most often they are asymmetrical short-circuits. Symmetrical short-circuits rarely occur in practice as compared with asymmetrical ones (only a few percent of the total number of cases). Over half of short-circuits in overhead high voltage lines are transient short circuits [1]. Asymmetrical PS operating conditions cause many unfavorable phenomena, i.a. in synchronous generators. It causes the necessity to limit the duration of such operating conditions [2, 3, 4].

Due to the difficulty in modeling asymmetrical operating conditions, symmetrical short-circuits are mainly analyzed in simulation investigations. Most specialized programs for the analysis of PS transient states allow for simulation of only symmetrical operating conditions. Therefore, it is reasonable to conduct research aimed at simulating the disturbance waveforms of selected quantities for different asymmetries occurring in PS. For this purpose, commonly used models of generating unit elements can be used, but some modifications should be introduced in them.

The analysis of asymmetrical PS operating states, including short-circuits, may help e.g. in better selection of power protection settings [5, 6]. The values and waveforms of different PS quantities may vary considerably depending on the type of the occurring asymmetry. Only effectively functioning protections can reduce the negative effects of disturbances, which allows reducing the consequences of faults occurring in PS.

The aim of the paper is a comparison and a harmonic analysis of the disturbance waveforms of selected quantities for different short-circuits with earth and clear to earth in a single-machine PS, consisting of a generating unit, a high-voltage transmission line and a bus. There was also analyzed the influence of the distance between the short-circuit location and the influence of including or neglecting effects of particular elements of the generating unit model.

Model of the analyzed PS

As a part of the conducted research, a mathematical model of the PS was developed in the Matlab-Simulink environment. In this model, using Configurable Subsystem blocks, it is possible to conveniently configure a specific model of the generating unit by selecting models of its individual elements.

In the carried out calculations, there were used: the GENROU [7, 8] synchronous generator model with subtransient asymmetry X”_d ≠ X”_q [3, 6, 9, 10] as well as the models of the static excitation system operating in the Polish Power System [7], the IEEEG1 steam turbine [7, 8] and the PSS3B system stabilizer [7].

When analyzing asymmetrical operating conditions of the PS, it is convenient to express the equations of stator currents and voltages, transmission line and bus with use of phase quantities. The Park transformation is used to relate the quantities in the d, q, 0 coordinate system to those in the phase A, B, C coordinate system [2, 3, 6, 7, 9, 10, 11, 12].

In the investigations, there were analyzed various asymmetries occurring in the transmission line. Fig. 1 shows a diagram of the analyzed PS for short-circuits with earth and clear to earth. On its basis, the model of the power line and the bus was developed.

Fig.1. Diagram of the analyzed PS for short-circuits with earth (a) and clear of earth (b)

Symbols in the figure: i_j – generator stator currents (phase quantities, j = A, B, C), v_j – generator stator voltages, v_bj– bus voltages, Z_j– complex impedances of the transmission line, Z_sj, Z_bj – impedance of the j-th phase line fragments during the short-circuit, I_fd – generator field current, v_d – voltage between the neutral points of the generator and the bus, i_bj – bus currents, v_s – voltage at the short-circuit location, W_G – generator neutral point ground switch.

This schematic diagram for short-circuits with earth (Fig. 1a) applies to both normal PS operation and a short circuit. In order to model a earth short-circuit in selected phases, one should assume zero voltages of the bus v_bj in these phases and change the impedances of the line proportionally:

where: l – relative distance of the short-circuit location in the transmission line from the generating unit, in relation to the length of the whole line.

In connection with the assumed omission of transformation voltages [2, 6, 12], on the basis of Fig. 1, there were determined the algebraic relations between currents and voltages in the generator stator and transmission line equations. For healthy phases:

for phases short-circuited to ground:

where: φ_j – phase angles of the particular transmission line impedances, f = 50 Hz. From (2a) and (2b) three equations were derived. In addition, for the system with ungrounded generator neutral point (open switch W_G in Fig. 1a):

Using the generator stator axial voltages (output signals of the generator model and the generating unit model) and the axial voltages of the bus, it is possible to determine the phase voltages by the inverse Park transformation. Hence, from the system of equations (2), one can calculate four unknows: three phase currents and possibly the voltage vd at successive time instants. This is implemented in the developed model in the Matlab-Simulink environment. Using the formulas for the generator phase currents and Park transformation, the generator axial currents, being the input signals of the generator model and the generating unit model, were determined.

Fig. 1b presents a diagram for a 2-phase short-circuit clear of earth. The following equations determine algebraic relations between currents and voltages for this fault:

where: Δt_j, Δt_bj– time delays of the current waveforms, defined as in (2a); other symbols as in formulas (2).

From the system of equations (3) one can calculate seven quantities: five phase currents and voltages: v_s and v_d. The thus determined PS model is complete and allows for carrying out simulation calculations.

Exemplary calculations

In the calculations presented, the generator worked with grounded neutral point (closed switch W_G in Fig. 1a, voltage v_d = 0) in order to check an effect of the zero axial components of stator currents and voltages on the phase waveforms. In real high-voltage power systems, generators work usually with isolated neutral points.

In the first case, simulation calculations were performed for 3-phase short-circuits and 2-phase short-circuits with earth (in phases A and B) with a duration time equal to 0.15 s. There were analyzed short-circuits at different distances l from the generating unit for the generating unit with and without a PSS. It was assumed that in the steady state before the short-circuit, the generator was loaded with active power P₀ = 0.5 p.u. and reactive power Q₀ = 0.2 p.u.

Figures 2 and 3 show the waveforms (in relative units p.u. [4, 6, 7, 10, 11, 12]) of the angular speed deviation of the generator rotor Δω and the generator instantaneous power P in the analyzed cases, for different relative distance l from the generating unit.

In the second case, there were analyzed long-lasting short-circuits: 1-phase with earth (in phase A), 2-phase clear of earth (in phases A and B) and 3-phase, at the distance l = 1% from the generating unit. In this case also breaks in the non short-circuited phases were assumed. In the steady state before the short-circuit the following load was assumed: P₀ = 0.1 p.u. and Q₀ = 0.05 p.u. The influence of the excitation system, the turbine and the PSS was neglected in this case.

Fig. 4 show the waveforms of stator current i_A (in relative units) for the analyzed fault types.

Tab. 1 presents the harmonic amplitude distributions of the current in phase A, the voltage in phase C and the field current in the steady state of the short-circuit. The higher harmonics percentage values are given in relation to the first harmonic for the stator quantities, and in relation to the constant component for the generator field current. The reference values in relative units (p.u.) are given in brackets.

Conclusions

The calculations made allowed formulating the following conclusions:

• The use of the transmission line and bus model for the phase components of currents and voltages in the developed PS model allows for easy modeling of various symmetrical and asymmetrical short-circuits with earth and clear of earth, as well as breaks in individual phases of the transmission line.

• The disturbance waveforms of the analyzed quantities in the case of a two-phase short circuit with earth differ significantly from those for a three-phase short-circuit. During the two-phase short-circuit, the waveforms of the instantaneous power of the generator have components with relatively high frequencies (in comparison with the frequency of electromechanical swings). Such components do not occur under symmetrical operating conditions, including the three-phase short-circuit. After the short circuit, in each of the analyzed cases, the generating unit returns to symmetrical operation and only low-frequency components associated with electromechanical phenomena occur. After the three-phase short-circuit, the amplitudes of the deviations from the steady values of the waveforms, especially of the angular speed, are much larger than those after the two-phase short-circuit.

Fig.2. Waveforms of the generator rotor angular speed deviation Δω for the system without the PSS: for two-phase short-circuit with earth (a), for three-phase short-circuit (b); for the system with the PSS: for two-phase short-circuit with earth (c), for three-phase short-circuit (d)

Fig.3. Waveforms of the generator instantaneous power P for the system without the PSS: for two-phase short-circuit with earth (a), for three-phase short-circuit (b); for the system with the PSS: for two-phase short-circuit with earth (c), for three-phase short-circuit (d)

Table 1. Harmonic amplitudes of the analyzed quantities in the steady state

Fig.4. Currents in phase A: a) envelopes, b) enlargement of one period in the steady state

• The use of the PSS with appropriately selected settings allowed for significant increase in the damping of electromechanical swings in the PS. As a result, the angular stability of the PS improved significantly and the time needed to return to the steady state after the disturbance was shortened. In the analyzed cases, the PSS did not have a significant influence on the amplitudes of the deviations from the steady values of the waveforms during and immediately after the short-circuit.

• In the cases under consideration, the influence of the distance between the short-circuit location and the generating unit on the waveforms of the analyzed quantities is large only during the short-circuit duration. After the short-circuit, it is significant only in the case of the three-phase short-circuit for the system without the PSS, and in other cases it is negligible.

• The maximal amplitude of the generator long-lasting short-circuit current in the steady state occurs at the 1- phase short-circuit. The current amplitude in the steady state at the 2-phase short-circuit is higher than that at the 3- phase short-circuit. It complies with the synchronous machine theory [3, 13].

• In the cases of the asymmetrical short-circuits, odd higher harmonics with significant amplitudes occur in the steady state waveforms of the short-circuit current and the stator voltage on the non short-circuited phase. The generator field current in steady state includes a constant component and even higher harmonics. The harmonics distributions in the short-circuit current and the non short-circuited phase voltage are similar for both types of asymmetrical short-circuits. Only the third harmonic of the voltage at the 1-phase short-circuit has a much lower amplitude. The harmonics distributions in the field current are similar for the both types of asymmetrical short-circuits.

• At the 3-phase symmetrical short-circuit, the stator current and voltage waveforms includes practically only the first harmonic, and in the field current waveform only the constant component occurs. The generator subtransient asymmetry does not cause the occurrence of higher harmonics in case of the machine symmetric operation. The presented power system model is also used in other investigations. They focus on various RL and XT-type synchronous generator models [7, 8] taking into account or neglecting the stator transformation voltages, with different input and output signals. A further extension of the model is planned.

REFERENCES

[1] Kacejko P., Machowski J, Zwarcia w systemach elektroenergetycznych [Short-circuits in power systems] (in Polish), WNT, Warszawa, 2009
[2] Concordia Ch., Synchronous Machines. Theory and Performance, John Wiley & Sons, Inc., New York, 1951
[3] Boldea I., Synchronous Generators, Taylor & Francis, 2015
[4] Pyrhonen J., Jokinen T., Hrabovcova V., Design of Rotating Electrical Machines, Wiley & Sons, Ltd, 2008
[5] Ungrad H., Winkler W., Wiszniewski A., Protection techniques in electrical energy systems, Mercel Dekker Inc., New York 1995
[6] Machowski J., Bialek J., Bumby J., Power System Dynamics. Stability and Control, John Wiley & Sons, Chichester-New York, 2008
[7] Paszek S., Berhausen S., Boboń A., Majka Ł., Nocoń A., Pasko M., Pruski P., Kraszewski T., Measurement estimation of dynamic parameters of synchronous generators and excitation systems working in the National Power System (in Polish), Monograph No. 504, Wydawnictwo Politechniki Śląskiej, Gliwice, 2013
[8] Boboń A., Paszek S., Pruski. P., Kraszewski T., Bojarska M., Computer-aided determining of parameters of generating unit models basing on measurement tests, Przegląd Elektrotechniczny, 5 (2011), pp. 17-21
[9] Berhausen S., Boboń A., Determination of high power synchronous generator subtransient reactances based on the waveforms for a steady state two-phase short-circuit, Applied Mathematics and Computation, 319 (2018), pp. 538–550
[10] Wang X., Song Y., Irving M., Modern Power Systems Analysis, Springer Boston, MA, 2008
[11] Krause P., Wasynczuk O., Sudhoff S., Pekarek S., Analysis of Electric Machinery and Drive Systems, third ed., Wiley-IEEE Press, 2013
[12] Gao J., Zhang L., Wang X., AC machine systems, Mathematical Model and Parameters, Analysis, and System Performance, Springer-Verlag, Berlin-Heidelberg, 2009
[13] Krause P.C., Analysis of electric machinery, McGraw-Hill, 1986

Authors: dr inż. Piotr Pruski, E-mail: Piotr.Pruski@polsl.pl, prof. dr hab. inż. Stefan Paszek, E-mail: Stefan.Paszek@polsl.pl, Politechnika Śląska, Wydział Elektryczny, Instytut Elektrotechniki i Informatyki, ul. Akademicka 10, 44-100 Gliwice.

Source & Publisher Item Identifier: PRZEGLĄD ELEKTROTECHNICZNY, ISSN 0033-2097, R. 96 NR 2/2020. doi:10.15199/48.2020.02.05

Arc Plasma Energy Evolvement in 60 kV Network Circuit Breakers

Published by Salah Belkhir¹, Abderrahmane Ziani¹, Hakim Azizi2, Hocine Moulai¹,
University of Science and Technology Houari Boumediene, Algiers, Algeria (1), University Ziane Achour, Djelfa, Algeria (2) ORCID: 1. 0000-0002-5106-2821

Abstract. The evolvement of the electric and energetic properties of electric arcs at the poles opening of a circuit breaker (CB) is described by nonlinear mathematical models. Most of these models are dimensional types that do not describe the interaction between the arc and the network during the interruption phase. This paper is aimed at the determination of the energy necessary for the arc creation at the opening of a high-voltage circuit breaker with a black box model. The advantage of this model consists of its ability to link the intrinsic characteristics of the arc to the extern blowing (quenching) power. Moreover, it provides fast and stable solving without needing for spatial dimensions of the breaker. The model is applied to a line circuit breaker for which experimental results are available in the literature. Two phases of arc quenching evolvement are evidenced: Constant energy phase followed by a decreasing energy one. The Kema-based model is found to be more accurate for online plasma quenching analysis and the obtained results agree well with experimental ones where the heat energy represents the dominating part.

Streszczenie. Ewolucję właściwości elektrycznych i energetycznych łuków elektrycznych przy otwarciu biegunów wyłącznika opisują nieliniowe modele matematyczne. Większość z tych modeli to modele wymiarowe, które nie opisują interakcji między łukiem a siecią podczas fazy przerwania. Celem artykułu jest określenie energii niezbędnej do wytworzenia łuku przy otwarciu wyłącznika wysokonapięciowego z modelem czarnej skrzynki. Zaletą tego modelu jest możliwość powiązania wewnętrznych charakterystyk łuku z zewnętrzną mocą nadmuchu (gaszenia). Ponadto zapewnia szybkie i stabilne rozwiązywanie bez konieczności wymiarowania przestrzennego wyłącznika. Model stosuje się do wyłącznika liniowego, którego wyniki eksperymentalne są dostępne w literaturze. Wykazano dwie fazy rozwoju gaszenia łuku: faza stałej energii, po której następuje faza malejącej energii. Stwierdzono, że model oparty na Kema jest dokładniejszy do analizy hartowania plazmowego w trybie online, a uzyskane wyniki dobrze zgadzają się z wynikami eksperymentalnymi, w których dominującą część stanowi energia cieplna. (Ocena energii łuku plazmowego w wyłącznikach sieciowych 60 kV)

Keywords: High voltage, circuit breaker, energy, plasma quenching.
Słowa kluczowe: Wysokie napięcie, wyłącznik, energia, hartowanie plazmowe.

Introduction

High voltage circuit breakers (CB) are necessary switchgear for the protection of electrical networks. In both low and high voltage networks, there is a very large number of breaking techniques that use the electric arc as a way to evacuate energy [1-5]. Also, the low cut-off times of less than some milliseconds and the high energy released during the electric arc formation in high voltage circuit breakers make the measurements difficult to achieve and also expensive [6-8]. The development and manufacture of HV circuit breakers require a detailed knowledge of heat transfer mechanisms that evolve at the arc extinction [2, 9, 10]. Although these phenomena are important and of great interest, one notes that only few articles are devoted to.

To extinguish the electric arc appearing at the opening of a high voltage circuit breaker, all the electromagnetic energy stored by the network must be dissipated [1, 11, 12]. Following the high Joule energy released during the CB opening which can reach 30000 J [11, 13, 14], the arc can be cooled through a strong blow of a suitable gas.

In this paper, a power balance has been set up in order to model the energy needed for the creation of the arc during the extinction phase. A precise description of that energy is essential because it determines the dielectric recovery of the arc. The model has been implemented in the Simulink- Matlab environment to determine at first the arc voltage and conductance, and then to calculate the arc creation energies. The model takes into account the ionization time constant of the arc, its voltage, current and even the external blowing power.

Mechanisms of the electric arc formation

Thanks to a mechanical system and as a result of fault current occurrence the high voltage circuit breaker electrodes separate in a quenching chamber. However, after contacts separation, an arc appears and the current continues to flow. In the presence of a gas, this arc is associated to the corresponding plasma.

The filling gas is usually SF6, chosen for its excellent thermal and dielectric properties [11]. The electric arc occurs in the zone of high ionic and electronic density provided from the inter-contacts medium or metal vapors from the circuit breaker poles [2]. The junction zones that bridge the arc column to the contacts are at temperatures close to the melting point of the metal, hence the thermoionic emission is possible [15]. Thus, the electric arc consists of plasma composed of ionized gas and metal vapors.

The arc quenching Models found in the literature [14, 16-20] are always based on Maxwell’s equations, Ohm’s law and the conservation equations of mass and energy. Delalondre et al. [21] developed a two-dimensional code to simulate switching in high voltage circuit breakers. They obtain fields of temperature and potential close to experimental ones. Chevrier et al proposed in [22] a switching arc model in low voltage circuit breakers, while Lindmayer et al [23] have developed a three dimensional model for a low-voltage circuit breaker which predicts the arc movement as a function of temperature and pressure. In addition, Gonzalez [14] adapts a commercial code (Fluent) for thermal plasma behavior in low-voltage circuit breakers. Cassie and Mayr, Lowke et al [24] and Wang et al. [25] have modeled two-dimensional variations of temperature and conductance of quenching arcs in air and SF6. On the other hand, Schavemaker et al [18] and Guardado et al [19] show that a 0D model is sufficient to follow the evolution of arc voltage in high voltage circuit breakers.

Energy balance of the electric arc

The rich bibliography about energy transfer available in [26-33] shows that energy models are based on the conservation equations of mass and energy.

where ρ and H are respectively the density and enthalpy of the plasma, v^→ its velocity, j^→ the current density, B^→ the magnetic field induction, P_ray the power lost by radiation and T the thermodynamic temperature.

Experimental measurements have shown that the term (j^→ Λ B^→) representing the induced current in the arc is often negligible [14], except for vacuum circuit breakers [34]. The appearance of the electric arc of current i(t) at the separation of the contacts creates an arc voltage U(t) that will determine a thermal elementary Joule energy dW during a time dt, thereby producing a very high temperature rise. The expression of this energy is given by:

The energy balance per time unit or the power balance reported to the arc in absence of magnetization is expressed by the following relationship:

Where Parc is the total electric power, P_J represents the Joule power provided to the arc which plays a role in the arc temperature rise, P is the cooling capacity due to the blowing of SF_6, P_C the power lost by thermal conduction and P_ray the power dissipated by radiation.

The term P_C can be expressed as a function of the arc temperature T and the temperature T₀ of the external environment by the following relationship:

where, K is the thermal conductivity of the arc and r its radius.

Even if their dissipation assessment is necessary [35, 36], the terms P_ray and P_C are usually neglected in black box models calculations [18]. Thus, P_ray and P_C are estimated to 1% in [21]. Rachard et al [29] reported that the thermal energy dissipated by conduction has moderate influence on the power balance.

New energy model approach

The arc models can be classified into two groups: physical models based on hydrodynamic codes deducted from Navier-Stokes and black box models where the local properties are averaged and their governing equations do not make appear differential terms on the space variables. This is why they are called 0D models. These include the models of Mayr, Cassie and Kema [11, 13].

Also, to follow the energy Q of arc creation in a circuit breaker during its extinction phase, the following assumptions will be adopted. They have been used by several authors [11, 13, 37] in 0D models, namely:

– The arc column is cylindrical in shape.

– The resistivity is constant and its section decreases during extinction.

– The electric field within the arc is constant – The energy of arc creation is proportional to its base surface.

To establish a simple model governing the arc formation, one must assume that the conductance g is only expressed as a function of energy Q used for this arc formation.

Q represents macroscopically the plasma ionization energy [11, 13].

Thanks to a differential equation, we can also write that the difference between the electric power supplied by the network u ⋅ i and the cooling power P injected by blowing is used to create the arc.

where P is the total cooling power provided to the arc, u the arc voltage, i the current through the arc and dt/dQ the power necessary for the arc creation. By expressing the differential of equation (5) as a function of time by multiplying and dividing by the same quantity dQ , we obtain:

Ohm’s law for a constant electric field E and a current i(t) crossing through an arc of cylindrical geometry and radius R gives for an electric resistivity ρ:

The conductance per unit length of a cylindrical arc can be expressed by:

where S is the surface of arc column. The arc section can be then deduced:

S= g x ρ

The surface of the arc assumed in the Cassie’s assumption [22] has been used to determine the energy Q necessary to create the arc:

One obtains then:

So, we have:

And by substituting the value of C_c by Q/S and the value of g by S/ρ, we obtain:

This equation can be replaced by the logarithmic expressions between Q and g:

Three 0D models (black box) are used to describe the evolution of the term Q through g. By using the Mayr equation, we obtain the following expression:

Where, P represents the blowing power in Watts and τ the deionization constant of the blowing gas (SF₆).

Similarly, Q can be followed by Cassie equation given by the following expression:

where u_c is the Cassie constant voltage. By using the Kema equation, one then obtains:

Where P is the external blowing power and P₁ a regulation coefficient that is 0.9943. The consecutive arc column surface evolvement is depicted on figure1.

Presentation of the used Electrical network

The HV electrical network used for simulation of breaking arc is shown in Figure 2. Such approach is suitable to take into account external parameters of the CB including the topology of the network and the connected loads [38, 39]. The modelled CB has the same characteristics than that used by Schavemaker et al. [18], namely longitudinal impedance including a resistance and reactance per unit length of the line and two transverse admittances. This network is powered by an electromotive force e = 60 kV. The network characteristics are as follows:

Inductance L=3.5×10^-3 H, Resistance R=30 Ω and Capacitance C=2.10^-6 F. The frequency is set to 50 Hz. Figure 2 shows the Matlab Simulink synoptic diagram of the studied network. The simulations were performed thanks to the SIMpower SYSTEM tool.

Fig. 1. Arc extinction physical behaviour

Simulation

In this work, we focused on a typical HV SF₆ circuit breaker used in substations and for which experimental results are available in the literature [13, 14, 16-19]. The simulation begins by initializing the arc parameters. The initial conductance g0 of the plasma is first fixed to10⁴ S/m which corresponds to a conducting state. It will allow us to compare our results with those obtained by Schavemaker [18]. The initial energy Q0 is then set to 10⁴ J.

Figure 2. Matlab Simulink synoptic diagram of the studied network

Figure 3. DEE block used to solve the Mayr Model

The circuit breaker opening time was set to 0.02 s. The blowing power of SF₆ was simulated by the Step block of Simulink with a step equal to the opening time of the circuit breaker. The simulation time was fixed between 0 and 0.03 s to allow comparisons with experimental results available in the literature. Several solvers of differential equations have been tested to obtain the best convergence of the solution to finally choose the ode45/Matlab solver with a variable time step in order to satisfy a relative tolerance of 10^-3.

Mayr, Cassie and Kema differential equations solving was made by the use of the Differential Equation Editor (DEE) block. The conductance was replaced by the variable x and the current by the expression u.exp(x). The solutions of the equations are generated in the DEE block as currents that are afterwards injected into the network. The steps followed to solve the Mayr model are presented below. Thus the differential equation was multiplied by the variable u(2) that has been introduced to control the opening time of the circuit breaker by providing a zero signal since the opening time is not reached.

By setting ln(g) = x , we obtain:

u(2)=0 for a time less than the breaker opening time.

u(2)=1 : For a time greater than or equal to the breaker opening time.

where g is the arc conductance, u(1) the arc voltage and i the arc current.

The output signal is a function of type y= exp(x(1)) x u(1) . To transform it into a usable current, a controlled current source is inserted at the output of the DEE. Figure 3 shows the adopted solving diagram in the Simulink environment.

Results and Interpretations

The variations of arc voltage in SF6 through the three models, namely Mayr, Cassie and Kema, for τ = 0.3μs have been plotted on figure 4. The results analysis shows that Mayr model is highly compatible with the Kema model, unlike Cassie model that reproduces constant arc voltages. One observes a voltage peak of about 80 kV which corresponds to the juxtaposition of the transient recovery voltage (TRV) with the numerically calculated voltage presenting good agreement with the experimental values measured by Schavemeker et al. [18] and Guardado et al. [19]. These authors have performed arc voltage measurements at the opening of a high voltage circuit breaker and also observed a sudden voltage increase at the opening of the poles. In Figure 5, the energies of arc creation during the extinction are presented through the model of Mayr. One can note a sudden decrease in energy during the first milliseconds following the arc creation. The shape of the obtained curves depends greatly on the external cooling power P where a fast decrease is observed for P = 30900 W. The numerical solution of Q based on Kema model for three powers is presented on Figure 6. This simulation shows two phases in the evolution of the arc. During the first phase, the creation energy remains constant, while in the second phase, a decrease of this energy is observed. One notes that the energy coupling with the Kema model is better adapted than with Mayr model. In addition, from the results reported in figure 7, i appears that the Cassie model is also not suitable for the simulation of Q. Figure 8 shows the joule heat energy injected during the opening of the circuit breaker poles. There is a rise of the energy necessary to maintain the arc. The injected energy is about 5×10⁵J for Kema and 1.4×10⁵J for Mayr model.

Figure 4. Arc voltage variations as a function of time for the considered models

Figure 5. Energy evolvement due to the arc versus time for the Mayr model

Figure 6. Energy evolvement due to the arc versus time for the Kema model

Figure 7. Energy evolvement due to the arc versus time for the Cassie model

By analysing the obtained results, it appears that the Q energy is about 1% higher than the values of the injected Joule thermal energy. These observations agree well with the results reported by Rachard et al [29] and Van Der Sluis et al. [31]. These authors show that the heat energy of the arc is dominating in high voltage circuit breakers.

Figure 8. Injected heat energy versus time obtained with Kema and Mayr models

Conclusion

A new model of arc creation energy at the opening of high-voltage circuit breakers was developed thanks to acceptable assumptions. Three conductance models were implemented and simulated in Matlab-Simulink environment to retain the best coupling between energy and conductance. The simulations have shown a decrease of the creation energy during the arc extinction phase.

The developed hybrid method has the advantage to no need for spatial dimensions to solve accurately the energy equations. Moreover, the equations are easy to implement and their solving is faster and more stable than the all numerical methods provided by commercial software.

The Kema-based model is found to be more accurate for online plasma quenching analysis where two phases of arc quenching evolvement are evidenced: Constant energy phase followed by a decreasing energy one. The obtained results are found to agree well with experimental ones where the heat energy represents the dominating part.

REFERENCES

[1] T. Chmielewski, P. Oramus, M. Szewczyk, T. Kuczek, W. Piasecki, Circuit breaker models for simulations of short-circuit current breaking and slow-front overvoltages in HV systems, Electric Power Systems Research, 143 (2017) 174-181. DOI: 10.1016/j.epsr.2016.10.046
[2] V. Abbasi, A. Gholami, K. Niayesh, The Effects of SF6-Cu Mixture on the Arc Characteristics in a Medium Voltage Puffer Gas Circuit Breaker due to Variation of Thermodynamic Properties and Transport Coefficients, Plasma Science and Technology, 15 (2013) 586-592. https://doi.org/10.1088/1009-0630/15/6/18
[3] F. Yang, Y. Wu, M.Z. Rong, H. Sun, A.B. Murphy, Z. Ren, C. Niu, Low-voltage circuit breaker arcs—simulation and measurements J. Phys. D: Appl. Phys. 46 (2013) 273001. https://doi.org/10.1088/0022-3727/46/27/273001
[4] M. Mürmann, A. Chusov, R. Fuchs, A. Nefedov, H. Nordborg, Modeling and simulation of the current quenching behavior of a line lightning protection device, J. Phys. D: Appl. Phys. 50 (2017) 105203. doi:10.1088/1361-6463/aa560e
[5] J. Valenta, M. Samohejl, M. Fendrych, P. Kloc, L. Dostál, Diagnostics of Various Phenomena in LV Devices under Real Switching Conditions, Plasma Physics and Technology 4 (2017) 257–260. doi:10.14311/ppt.2017.3.257
[6] T. Cheng, W. Gao, W. Liu and R. Li, Evaluation method of contact erosion for high voltage SF6 circuit breakers using dynamic contact resistance measurement, Electric Power Systems Research, 163 (2017) 725-732. DOI: 10.1016/j.epsr.2017.08.030
[7] D. Dufournet, G.F. Montillet, Transient recovery voltages requirements for system source fault interrupting by small generator circuit breakers, IEEE Trans. Power Delivery, 17 (2002) 474-478. DOI: 10.1109/61.997921
[8] E.A.L. Vianna, A.R. Abaide, L.N. Canha, V. Miranda, Substations SF6 circuit breakers: Reliability evaluation based on equipment condition, Electric Power Systems Research, 142 (2017) 36-46. DOI: 10.1016/j.epsr.2016.08.018
[9] Y. Wu, H. Sun, Y. Tanaka, K. Tomita, M. Rong, F. Yang, Y. Uesugi, T. Ishijima, X. Wang, Y. Feng, Influence of the gas flow rate on the nonchemical equilibrium N2 arc behavior in a model nozzle circuit breaker, J. Phys. D: Appl. Phys., 49 (2016) 425202. https://doi.org/10.1088/0022-3727/49/42/425202
[10] Y. Wu, Y. Cui, J. Duan, H. Sun, C. Wang and C. Niu, Influence of arc current and pressure on non-chemical equilibrium air arc behavior, Plasma Science and Technology, 20 (2017) 014021. https://doi.org/10.1088/2058-6272/aa9325
[11] S. Vacquié, Arc électrique, Techniques de l’Ingénieur, Traité de Génie Electrique D 2870, 1995.
[12] V. Abbasi, A. Gholami, K. Niayesh, Impact of radial external magnetic field on plasma deformation during contact opening in SF6 circuit breakers, Journal of Physics D: Applied Physics, 45 (2012) 415201. https://doi.org/10.1088/0022-3727/45/41/415201
[13] CIGRE Working Group 13-01, Applications of Black Box Modelling to Circuit Breaker, Electra, 149 (1993) 41-71.
[14] J.J. Gonzalez, F. Lago, P. Freton, M. Masquère, X. Franceries, Numerical modelling of an electric arc and its interaction with the anode : Part ii. the three-dimensional model-influence of external forces on the arc column., J. Phys. D : Appl. Phys., 38 (2005) 306–318. https://doi.org/10.1088/0022-3727/38/2/016
[15] X. Liu, S. Wang, Y. Zhou, Z. Wu, K . Xie, N. Wang, Thermal radiation properties of PTFE plasma, Plasma Science and Technology, 19 (2017) 064012. https://doi.org/10.1088/2058-6272/aa65e8
[16] A. Ziani, H. Moulai, Extinction properties of electric arcs in high voltage circuit breakers, Journal of Physics D: Applied Physics, 42 (2009) 105205. https://doi.org/10.1088/0022-3727/42/10/105205
[17] A. Ziani, H. Moulai, 0D Model of thermal exchanges at the opening of an SF6 high voltage circuit breaker, ACTA Press, Proc. of the Int. Conf. on Power and Energy Systems, September 7–9, 2009, Palma de Mallorca, Spain, paper 681-056. http://www.actapress.com/Abstract.aspx?paperId=35434
[18] P.H. Schavemeker, L. Van der Sluis, An improved Mayr type arc model based on current zero measurement, IEEE Trans. Power Delivery, 15 (2000) 580-584. DOI: 10.1109/61.852988
[19] J.L. Guardado, S.G. Maximov, E. Melgoza, J.L. Naredo, P. Moreno, An Improved Arc Model Before Current Zero Based on the Combined Mayr and Cassie Arc Models, IEEE Trans. Power Delivery, 20 (2005) 138-142. DOI: 10.1109/TPWRD.2004.837814
[20] J.B. Belhaouari, Modélisation de l’extinction d’un arc de SF6 hors d’équilibre thermodynamique local, Doctorat thesis, Paul Sabatier university, Toulouse III, France, N°2780, 1997. tel-00003150v2
[21] O. Simonin, C. Delalondre, P.L. Viollet, Modelling in thermal plasma and electric arc column, Pure and Appl. Chem., 64 (1992) 623–628,. DOI: 10.1351/pac199264050623
[22] P. Chévrier, M. Barrault, C. Fiévet, J. Maftoul, J.M. Frémillon, Industrial applications of high-, medium- and low-voltage arc modelling. J. Phys. D: Appl. Phys., 30 (1997) 1346–1355. https://doi.org/10.1088/0022-3727/30/9/010
[23] M. Lindmayer, E. Marzahn, A. Mutzke, M. Springstubbe, Lowvoltage switching arcs – experiments and modeling, 15th Symposium on Physics of Switching Arc, Brno, Czech Republic, 2003.
[24] J.J. Lowke, R.E. Voshall, H.C. Ludwing, Decay of electrical conductance and temperature of arc plasmas, J Appl. Physics, 44 (1973) 3513–3523. DOI: 10.1063/1.1662795
[25] W.Z. Wang, J.D. Yan, M.Z. Rong, A.B. Murphy, J.W. Spencer, Theoretical investigation of the decay of an SF6 gas-blast arc using a two-temperature hydrodynamic model, Journal of Physics D: Applied Physics, 46 (2013) 065203. https://doi.org/10.1088/0022-3727/46/6/065203
[26] M. Razafinimanana, A. Gleizes, F. Mbolidi, S. Vacquié, D. Gravelle, Experimental study of an SF6 arc in extinction, J. Phys. D: Appl. Phys, 23 (1990) 1671. https://doi.org/10.1088/0022-3727/23/12/027
[27] J.C. Lee, Y.J. Kim, SF6 arc plasma modelling for compact and environmental-friendly gas circuit breaker, Surface and Coatings Technology, 201 (2007) 5641-5645. DOI: 10.1016/j.surfcoat.2006.07.110
[28] A. Ziani, H. Moulai, Thermal Transfers of SF6 Electrical Arcs in High Voltage Circuit Breakers, Acta Physica Polonica A, 123 (2012) 241-244. DOI: 10.12693/APhysPolA.123.241
[29] H. Rachard, P. Chévrier, D. Henry, D. Jeandel, Numerical study of coupled electromagnetic and aerothermodynamic phenomena in circuit breaker electric arc, Int. J. Heat and Mass Transfer, 42 (1999) 1723-1734. DOI: 10.1016/S0017-9310(98)00110-0
[30] A. Gleizes, A.M. Rahal, H. DeLacroix, P. Van Doan, Study of a circuit-breaker arc with self-generated flow. I. Energy transfer in the high-current phase, IEEE Transactions on Plasma Science, 16 (1989) 606 – 614. DOI: 10.1109/27.16548
[31] L.Van Der Sluis, W.R. Rutgers, C.G.A. Koreman, A Physical Arc Model for the Simulation of Current Zero Behaviour of High-Voltage Circuit Breakers, IEEE Trans. on Power Delivery, 7 (1992) 1016-1022. DOI: 10.1109/61.127112
[32] N. Osawa, Y. Yoshioka, Analysis of nozzle ablation characteristics of gas circuit breaker, IEEE Trans. Power Del., 25 (2003) 810 – 815. DOI: 10.1109/TDC.2003.1335379
[33] I.M. Dudurych, T. J. Gallagher, E. Rosolowski, Arc effect on single-phase reclosing time of a UHV power transmission line, IEEE Trans. Power Delivery, 19 (2004) 854-860. DOI: 10.1109/TPWRD.2004.824404
[34] Y. Chen, F. Yang, H. Sun, Y. Wu, C. Niu, M. Rong, Influence of the axial magnetic field on sheath development after current zero in a vacuum circuit breaker, Plasma Science and Technology, 19 (2017) 064003. DOI: 10.1088/2058-6272/aa65c8
[35] C. Jan, Y. Cressault, A. Gleizes, K. Bousoltane, Calculation of radiative properties of SF6–C2F4 thermal plasmas—application to radiative transfer in high-voltage circuit breakers modeling, Journal of Physics D: Applied Physics, 47 (2014) 5204. DOI: 10.1088/0022-3727/47/1/015204
[36] S. Tsuda, K. Horinouchi, H. Yugami, Enhancing the radiative heat dissipation from high-temperature SF6 gas plasma by using selective absorbers, Journal of Physics D: Applied Physics, 50 (2017) 365601. DOI: 10.1088/1361-6463/aa7fd5
[37] A. Ziani, H. Moulai, Hybrid model of electric arcs in high voltage circuit breakers, Electric Power Systems Research, 92 (2012) 37-42. DOI: 10.1016/j.epsr.2012.04.021
[38] Joanna Budzisz, The model of a vacuum circuit breaker for switching on capacitor bank, Przegląd Elektrotechniczny, ISSN 0033-2097, R. 95 NR 2/2019, doi:10.15199/48.2019.02.31.
[39] Joanna Budzisz, Zbigniew Wróbleski, The model of a vacuum circuit breaker in MATLAB software for the analysis of overvoltages and overcurrents in capacitive electrical circuits, Przegląd Elektrotechniczny, ISSN 0033-2097, R. 92 NR2/2016, doi:10.15199/48.2016.02.37

Authors: Dr. Salah Belkhir, University of Science and technology Houari Boumediene, Bab Ezzouar, Algiers 16025 Algeria; Email: belkhir.s@hotmail.com; Prof. Abderrahmane Ziani, University of Science and technology Houari Boumediene, Bab Ezzouar, Algiers 16025 Algeria; Email: ziani08@yahoo.fr; Dr. Hakim Azizi, University Ziane Achour, Djelfa, 17000 Algeria; Email: azizihakimbilal@yahoo.fr; Prof. Hocine Moulai, University of Science and technology Houari Boumediene, Bab Ezzouar, Algiers 16025 Algeria; Email: hmoulai@usthb.dz

Source & Publisher Item Identifier: PRZEGLĄD ELEKTROTECHNICZNY, ISSN 0033-2097, R. 98 NR 10/2022. doi:10.15199/48.2022.10.09

Comparison of the Results of Simulation Modeling of an Asynchronous Electric Motor with the Calculated Electrodynamic and Energy Characteristics

Published by 1. Oleg GUBAREVYCH¹, 2. Svitlana GOLUBIEVA¹, 3. Inna MELKONOVA², State University of Infrastructure and Technologies (1), Volodymyr Dahl East Ukrainian National University (2)
ORCID: 1. 0000-0001-7864-0831; 2. 0000-0001-8285-7566; 3. 0000-0001-6173-1470

Abstract. The paper presents the results of a comparison of the electrodynamic and energy characteristics and parameters of an asynchronous motor, obtained by simulation and calculated by the classical method. The mathematical model in the MATLab software environment is used for research. The research results are relevant when choosing and using the proposed simulation model of three-phase squirrel-cage asynchronous motors for further research, including the effect of various engine defects on its performance

Streszczenie. W artykule przedstawiono wyniki porównania charakterystyk i parametrów elektrodynamicznych i energetycznych silnika asynchronicznego, uzyskanych metodą symulacji i obliczonych metodą klasyczną. Do badań wykorzystano model matematyczny wykonany w środowisku oprogramowania MATLab. Wyniki prac mają istotne znaczenie dla wyboru i wykorzystania zaproponowanego modelu symulacyjnego asynchronicznego silnika elektrycznego z wirnikiem klatkowym do dalszych badań, w tym wpływu różnego rodzaju uszkodzeń silnika na jego pracę. (Porównanie wyników modelowania symulacyjnego asynchronicznego silnika elektrycznego z obliczonymi charakterystykami elektrodynamicznymi i energetycznymi)

Keywords: asynchronous motor, simulation modeling, charakterystyka elektrodynamiczna, mathematical model.
Słowa kluczowe: silnik asynchroniczny, modelowanie symulacyjne, electrodynamic characteristics, model matematyczny

Introduction

Improving the operational reliability of electromechanical equipment is a modern priority task, the solution of which determines the result of the efficient operation of industrial and transport enterprises. Three-phase squirrel-cage asynchronous motors are among the most common electrical machines used to drive various mechanisms in all industries. Currently, almost 70% of the machines used in industry are three-phase asynchronous motors because they are simple, reliable and inexpensive [1]. The requirements for each technological process determine the need to set and maintain the operating parameters of the motors used with high accuracy at a given level. Many malfunctions and defects that occur in difficult operating conditions quickly progress and disable electric motors, even with a short service life, leading them to an emergency stop. Timely and reliable detection of damage not only increases the reliability of motors, but significantly reduces repair time and reduces unforeseen costs. Experience in the operation of asynchronous electric motors shows that the development and implementation of modern diagnostic tools based on a thorough study of the processes occurring when defects occur is one of the most important and effective factors in increasing the economic efficiency of using electromechanical equipment in industry [2, 3].

Determining the type and degree of damage, establishing their influence on the performance of electric motors, improving the accuracy of predicting the final resource and uptime is possible on the basis of studying electrodynamic processes that occur in the presence of defects of various kinds [4-6]. Mathematical modeling methods are widely used to conduct research in various fields [7-9]. Simulation models of asynchronous electric motors for the correct use of the results in improving diagnostic systems and studying processes occurring at different degrees of defects must be checked for compliance with the properties of the simulated object to real processes. The results of simulation modeling are of an estimated nature and require further verification, in particular, by conducting full-fledged experimental studies on real objects. Given that experimental studies are highly laborious, the paper proposes an approach to evaluating the results of simulation modeling by comparing them with a set of calculated energy and electrodynamic indicators of an asynchronous motor.

The aim of the article is to conduct a study on the evaluation of the selected mathematical model of an asynchronous motor by com-paring the results of simulation modeling and the electrodynamic characteristics and energy parameters calculated by the classical method. The use of a mathematical model with the established accuracy of the results of simulation modeling in further research will allow more correct consideration of the influence of defects in the operation of asynchronous motors to determine the method for their diagnosis and assessment of the degree of damage and study of the ongoing electrodynamic processes.

To achieve the aim, the following tasks were completed:

– the parameters and performance characteristics of the selected base motor were calculated according to the classical method;

– simulation modeling in the MATLab software environment of the operation of an asynchronous motor was carried out using the selected mathematical model;

– obtained as a result of mathematical modeling electrodynamic characteristics and energy parameters of an asynchronous electric motor;

– a comparison of the results of modeling and calculations using the electrodynamic characteristics of an asynchronous electric motor was carried out.

Choice of a model for the study of electrodynamic processes and the basic type of motor

The statistics of operating experience of asynchronous motors shows that the largest share of operability failures is due to failures in the stator and, according to various data, taking into account the operating area, 77-85%, the rotor accounts for 6-8% and the bearing assembly – up to 8-14% [1, 3, 10]. The main part of the defects in the motor leads to the occurrence of an asymmetric rotating stator field. Therefore, to conduct a study on the manifestation of a larger range of damage, it is necessary to use a reliable mathematical model of an asynchronous motor with the possibility of research, including under asymmetric modes that occur during operation with various types of stator defects [11-13].

To conduct research on electromagnetic processes, there are a large number of approaches to simulation modeling of asynchronous motors. Their differences are mainly related to the choice of the coordinate system in which differential equations are composed that describe the operation of an asynchronous motor. In works [6, 11], a model is considered in which the equations de-scribing the operation of an asynchronous motor are written in d-q coordinates, i.e. in a single-phase coordinate system. When using such a model, it becomes difficult to determine some parameters, for example, the imbalance of phase currents, which is necessary for diagnosing a number of defects. When modeling asynchronous motors with asymmetric windings, it is advisable to use a system of differential equations in hindered coordinates, as noted in [14]. To solve the problem, the use of other coordinate systems is incorrect. This is evidenced by studies in [15, 16].

To implement the simulation modeling of the operation mode of an asynchronous motor with asymmetrical windings, which occurs in the event of damage to one or more stator windings, one should set a change in the leakage inductance and active resistance of the corresponding winding (windings). That is, it is necessary to establish how much the actual values of the specified parameters differ from the nominal values. After that, take into account how the mutual inductance of the windings will change. To determine the change in the mutual inductance of the windings, it is necessary to establish what effect the change in the complex resistance of one winding (several windings) has on the inductance of the magnetic circuit. The papers [17, 18] show the established relationship between the winding inductances and the geometric dimensions of the windings, which must be taken into account when simulation modeling.

Thus, in order to conduct further research with a wider range of possible defects that affect the operating modes of motors, it is necessary to use a mathematical model of an asynchronous motor with the possibility of creating an asymmetric rotating field, made in “braking coordinates”, taking into account losses in steel and mechanical losses. In addition, for the implementation of this modeling principle, one should take into account the mutual inductance of the windings when the complex resistance of one or more phases of the stator winding, simulating its damage, changes. When determining the change in the mutual inductances of the windings, it is necessary to determine what effect the change in the complex resistance of one winding (two windings) has on the inductance of the magnetic circuit and establish the relationship between the inductances of the windings and the geometric dimensions of the windings, as well as the effect on the leakage inductance of each phase and mutual phase inductances.

The simulation model of an asynchronous motor proposed and used in the work, for which the reliability of real processes is established, is made in “braked coordinates”. The general view of the model and its implementation in the MATLab software environment are presented and discussed in detail in [19]. The implementation of the mathematical model [19] based on the configuration of mutual inductance with a change in the complex resistance of one or more phases of the motor is considered in [20]. This model can be adapted to study the operation of an asynchronous motor when such a defect occurs in it as an interturn short circuit of the stator windings, which entails an asymmetric rotation field, as well as to determine the starting and operating characteristics of the motor, calculate energy indicators when the asynchronous motor is operating with the specified defect. Thus, simulation modeling of an asynchronous electric motor was carried out using a mathematical model given and discussed in detail in [20].

The simulation model of an asynchronous motor, made in the MATLab software environment, is shown in fig. 1. The following parameters are displayed: stator phase voltages, stator phase currents, rotor phase currents, motor shaft speed and useful motor shaft torque. For this purpose, an oscilloscope implemented on the Scope element was used. This oscilloscope has four sections: U_s, I_s, I_r, n, M. The first section displays the voltages of the stator phases, the second – the phase currents of the stator, the third – the phase currents of the rotor, the fourth – the frequency of rotation of the motor shaft and the torque on the motor shaft. The signals corresponding to the torque on the motor shaft and the frequency of rotation of the motor shaft are displayed on the display of the Display measuring unit. This is necessary to determine the exact value of the indicated values. If it is necessary to measure the amplitudes and phase angles of the stator voltages, rotor currents, the Complex to Magnitude-Angle unit is used, at the input of which a response signal is applied, at the output there are the signals corresponding to the amplitude and phase of this signal. After that, the received signals are displayed on an indication organized using Display units.

Fig.1. Simulation model of an asynchronous motor

An asynchronous motor with a squirrel-cage rotor of the AIR132M4 type with a power of 11kW with a supply voltage of 220/380 V was chosen as the base motor for research. The rating data of the motor are given in Table 1.

For the nominal mode of the motor (see Table 1) according to the method given in [21], the following parameters were calculated:

– moment on the motor shaft; frequency of rotation of the motor shaft; useful power;

Table 1. The parameters of the sensor

– active, reactive and apparent power consumed from the network; losses in steel and copper of the stator and losses in the rotor; mechanical losses;

– phase current of the stator winding; Efficiency and power factor cosφ₁.

The supports are also calculated: the active resistance of the stator winding and the active resistance of the rotor winding, reduced to the stator winding; reactive supports of the stator winding and the resistance of the rotor winding, reduced to the stator winding and in the magnetizing circuit. Some calculated parameters of the prototype motor according to the classical method differ from the passport data. So, the error in calculating the useful power was 0,045%, the error in calculating the efficiency – 0,114%, the error in calculating cosφ₁ – 0,237%. The given deviations have the maximum differences among other calculated parameters, which is explained by the error of the calculation method used.

To evaluate the results of simulation modeling on the selected mathematical model of an asynchronous electric motor, according to the specified classical technique [21], the performance characteristics are calculated, shown in Table 2. The table highlights the values of the parameters corresponding to the nominal mode.

Results of simulation modeling of electrodynamic actions of an asynchronous electric motor

When setting symmetrical phase voltages of the stator on the model, the timing diagrams of which are shown in fig. 2, the values of the stator phase currents (fig. 3), the rotor phase currents (fig. 4) are obtained. The values of the stator phase voltage and the stator and rotor phase currents are given by the instantaneous values of these parameters.

Table 2. Calculation of the performance characteristics of an asynchronous motor

Fig.2. Timing diagrams of stator phase voltages for nominal mode

Fig.3. Timing diagrams of stator phase currents for nominal mode

Fig.4. Timing diagrams of rotor phase currents for nominal mode

To evaluate the results of the selected mathematical model, simulation modeling of the electrodynamic processes of the operation of an asynchronous electric motor was carried out with the passport data given in Table 1. The results of modeling the dependence of active power (P₁), useful power (P₂) consumed from the network, average stator current (I_1mid), efficiency (η), power factor (cosφ) and torque on the motor shaft (T) on the rotor speed (n₂) are given in Table 3. The simulation modeling was carried out with a fixed speed of the motor shaft.

Table 3. Results of simulation modeling of an asynchronous electric motor operation

Useful power in table 3 is determined by:

where T – the moment on the motor shaft, N·m, ω – the rotation speed of the motor shaft.

The rotation speed of the motor shaft for pairs of field-owls p=2 is determined by:

where n – the actual speed of the motor shaft, rpm.

The instantaneous reactive power is determined by:

Instantaneous apparent power consumed from the network:

Then the value of the instantaneous active power consumed from the network can be determined as:

The efficiency is determined by:

The power factor can be determined using the expression:

According to the results of simulation modeling (mod) given in Table 3 and calculated data (calc) from Table 2, the mechanical characteristics n₂=f(T) and the dependence of the rotor speed on the useful load on the shaft n₂=f(P₂) are constructed in fig. 5 and fig. 6, respectively.

Fig.6. Dependence of the rotor speed on the useful load on the shaft

Based on the results of tables 3 and 2, the dependences are plotted for the relative value of the active power consumed from the network on the relative value of the useful power P₁*=f(P₂*) (P₁/P_1н=f(Р₂/Р_2н) (Fig. 7) and the dependence of the average relative value phase current of the stator on the relative value of useful power (I_1mid*)=f(P₂*) (I₁/I_1н)=f(Р₂/Р_2н) (Fig. 8).

Fig.7. Dependence of the relative value of active power consumed from the network *(P₁*) (P₁*=P₁/P_1rat)* on the relative value of useful power (P₂*) (P2*=P₂/P_2rat)

Fig.8. Dependence of the average relative value of the stator phase current *(I_1mid*) (I_1mid*=I_1mid/I_1mid.rat)* on the relative value of the useful power *(P₂*) (P₂*=P₂/P_2rat)*

Based on the results of tables 3 and 2, the dependences of the energy indicators of the motor are plotted: efficiency η =f(Р₂/Р_2н) (Fig. 9) and power factor cosφ=f(Р₂/Р_2н), (Fig. 10) on the relative value of the useful power.

Fig.9. Dependence of the efficiency factor (η) on the relative value of the useful power *(P₂*) (P₂*=P₂/P_2rat)*

Fig.10. Dependence of the power factor (*cosφ*) on the relative value of the useful power *(P₂*) (P₂*=P₂/P_2rat)*

Analysis of the results of simulation modeling in comparison with the calculated electrodynamic and energy characteristics and motor parameters shown in figures 5-10 indicates a high degree of compliance of the results obtained with the calculated ones and a sufficiently high dynamic stability of the model in the working range (T=23,0- 110,0 N·m) (see fig. 5, 6). The greatest deviations are the discrepancies in the energy indicators of the motor, which are 1,2% for the efficiency factor and 0,8% for the power factor cosφ, which is explained by the errors of the calculation method.

Conclusion

Defects or damage to squirrel-cage asynchronous motors require timely diagnosis and study of their effect on the parameters and characteristics of operating electrical equipment. An effective means of assessing the influence of defects on the energy and electrodynamic processes occurring in an asynchronous motor is simulation modeling. Verification and evaluation of results of simulation modeling obtained on specific mathematical models for compliance with real ongoing processes is a necessary issue that helps to increase the level of correctness of the results when they are used in further studies of asynchronous electric motors. The most reliable results on establishing the level of model adequacy can be obtained by conducting experimental studies on real objects. Taking into account the significant laboriousness of conducting experimental studies in the work, an approach is proposed for evaluating the results of simulation modeling by comparing them with a set of calculated energy and electrodynamic indicators of an asynchronous motor, which is the novelty of this work.

Based on the results of the studies, it is found that the greatest discrepancies between the calculated data and the data obtained during modeling take place for a set of energy indicators and are: for efficiency factor – 1,2%, for power factor cosφ – 0,8%, which is explained by errors calculation method. Comparison of the electrodynamic characteristics and parameters shows a fairly high accuracy of the results obtained in the simulation modeling. Thus, the discrepancies in the torque values at the nominal rotor speed n₂ =1450 rpm are: the calculated value T=72,47 Nm, the value obtained by modeling T=72,443 Nm. The discrepancies in the values of the calculated and obtained by modeling torques over the entire operating range of the electric motor, with the rotation of the rotor shaft n₂=1417- 1485 rpm, which corresponds to T=23,0–110,0 Nm, also have an insignificant level that can be neglected.

The research results have shown the possibility of using the proposed simulation model for further research, in particular, the impact on the manifestation of various types of electrical defects in an asynchronous motor on its performance, with a fairly high correspondence of the results to the calculated values.

Subsequent work will be aimed at studying the effect of interturn short circuit in the phases of the stator winding of an asynchronous electric motor on its performance using the proposed mathematical model while improving the system for diagnosing asynchronous motors.

REFERENCES

[1] Khechekhouche A., Cherif H., Menacer A., Chehaidia S.E., Panchal H. Experimental diagnosis of inter-turns stator fault and unbalanced voltage supply in induction motor using MCSA and DWER. Periodicals of Engineering and Natural Sciences (PEN), 2020, vol. 8, no. 3, pp. 1202-1216
[2] Ciprian H., Szabó L. Wavelet Analysis and Park’s Vector Based Condition Monitoring of Induction Machines. / H. Ciprian, L. Szabó // Juornal of Computer Science and Control Systems, vol. 4, no. 2, 2011, pp. 35-38
[3] Gubarevych О.V, Goolak S.О., Golubieva S.M. An integrated approach to diagnosing asynchronous electric motors of water transport. New technologies, 2019, no. 2(9), pp. 48-61
[4] Mairte J., Gaboury S., Bbouchard B., Bouzouane, A. A new computational method for stator faults recognition in induction machines based on hyper-volumes. In: 2015 IEEE International Conference on Electro/Information Technology (EIT). IEEE, 2015, pp. 216-220
[5] Culbert I., Letal J. Signature analysis for online motor diagnostics: Early detection of rotating machine problems prior to failure, IEEE Industry applications magazine, 2017, no. 23(4), pp. 76-81
[6] Liubarskyi B., Petrenko O., Shaida V., Maslii A. Analysis of optimal operating modes of the induction traction drives for establishing a control algorithm over a semiconductor transducer. Eastern-European Journal of Enterprise Technologies, 2017, vol. 4, no. 8 (88), pp. 65-72
[7] Lovska A., Fomin O.A New Fastener To Ensure The Reliability Of A Passenger Car Body On A Train Ferry. Acta Polytechnica, 2020, vol. 60, no. 6, pp. 478–485
[8] Shavelkin A.A., Gerlici J., Shvedchykova І.О., Kravchenko K., Kruhliak H.V. Management of power consumption in a photovoltaic system with a storage battery connected to the network with multi-zone electricity pricing to supply the local facility own needs. Electrical Engineering and Electromechanics, 2021, no. 2, pp. 36–42
[9] Shavolkin O., Shvedchykova I., Demishonkova S., Pavlenko V. Przeglad. Increasing the efficiency of hybrid photoelectric system equipped with a storage battery to meet the needs of local object with generation of electricity into grid. Elektrotechniczny, 2021, no. 97(11), pp. 144–149
[10] Eldeeb H.H., Berzoy A., Saad A.A., Mohammed O.A. On-line Monitoring of Stator Inter-Turn Failures in DTC driven Asynchronous Motors using Mathematical Morphological Gradient. In 2019 IEEE Applied Power Electronics Conference and Exposition (APEC), 2019, pp. 1018-1023
[11] Ayyappan G.S., Ramesh Babu B., Srinivas K., Raja Raghavan M., Poonthalir R. Mathematical Modelling and IoT Enabled Instrumentation for Simulation & Emulation of Induction Motor Faults. IETE Journal of Research, 2021, pp. 1-13. doi: 10.1080/03772063.2021.1875272
[12] Ali M.Z., Shabbir M.N.S.K., Liang X., Zhang Y., Hu T. Machine Learning-Based Fault Diagnosis for Single-and Multi-Faults in Induction Motors Using Measured Stator Currents and Vibration Signals. IEEE Transactions on Industry Applications, 2019, no. 55(3), pp. 2378-2391
[13] Goolak S., Gerlici J., Gubarevych O., Lack T., Pustovetov M. Imitation Modeling of an Inter-Turn Short Circuit of an Asynchronous Motor Stator Winding for Diagnostics of Auxiliary Electric Drives of Transport Infrastructure. Communications – Scientific Letters of the University of Zilina, vol. 23 no. 2, pp. 65-74. doi: 10.26552/com.C.2021.2. pp. 65-74
[14] Pustovetov M, Soltus К., Sinyavskiy I. Computer simulation of induction motors and transformers. Examples of interaction with power electronic converters. LAP LAMBERT. Academic Publishing, 2013
[15] Yu M., Zhu J., Qiang D., Zhu Y. Numerical calculation of global temperature field during phase failure of small induction motor. In 2019 Chinese Control Conference (CCC), 2019, pp. 7143-7148
[16] Singh A., Grant B., DeFour R., Sharma C., Bahadoorsingh S. A review of induction motor fault modeling. Electric Power Systems Research, 2016, no. 133, pp. 191-197
[17] Goolak S. Methodical recommendations for application of the model of physical processes in a three-phase asynchronous motor. Proceedings of the State University of Infrastructure and Technology. Series: Transportation Systems and Technologies, 018, no. 1(32), pp. 4-13
[18] Goolak S., Gerlici J., Sapronova S., Tkachenko V., Lack T., Kravchenko K. Determination of Parameters of Asynchronous Electric Machines with Asymmetrical Windings of Electric Locomotives. Communications-Scientific letters of the University of Zilina, 2019, no. 21(2), pp. 24-31
[19] Pustovetov M.Yu. Approach to Computer Implementation of Mathematical Model of 3-Phase Induction Motor. IOP Conf. Series: Materials Science and Engineering [online]. 2018, no. 327(2), 022085
[20] Goolak S., Gubarevych О., Yermolenko E, Slobodyanyuk M, Gorobchenko O. Development of mathematical model of induction motor for vehicles. Eastern-European Journal of Enterprise Technologies. 2020, no. 2/2 (104), pp.24-35
[21] Kopylov I.P. Designing of electrical machines. Textbook for high schools. M.: Yurayt, 2019. 828р

Authors: PhD, Associate Professor Oleg Gubarevych, State University of Infrastructure and Technologies, st. Kyrylivska 9, 04071, Kyiv, Ukraine, E-mail: oleg.gbr@ukr.net; Svitlana Golubieva, State University of Infrastructure and Technologies, st. Kyrylivska 9, 04071, Kyiv, Ukraine, E-mail: glbvvnu@gmail.com; PhD, Associate Professor Inna Melkonova, Volodymyr Dahl East Ukrainian Nationa University, Prospect Tsentralnyiy 59a, 93400, Severodonetsk, Ukraine, E-mail: melkonova@snu.edu.ua,

Source & Publisher Item Identifier: PRZEGLĄD ELEKTROTECHNICZNY, ISSN 0033-2097, R. 98 NR 10/2022. doi:10.15199/48.2022.10.11

Determination of Electrical and Efficiency Parameters of Air Cooling of Low-Temperature PEM Fuel Cell Stack with Power of 5kW

Published by Andrzej RAŹNIAK¹, Magdalena DUDEK¹, Tomasz SIWEK¹, Piotr DUDEK², Wojciech KALAWA¹, AGH – Akademia Górniczo – Hutnicza im. Stanisława Staszica w Krakowie, Wydział Energetyki i Paliw (1) AGH – Akademia Górniczo – Hutnicza im. Stanisława Staszica w Krakowie, Wydział Inżynierii Mechanicznej i Robotyki (2)

Abstract: The possibilities of using low-temperature hydrogen-oxygen fuel cells with PEMFC proton-exchange membrane in transport and aviation were characterized and described. In this paper the emphasis was put on the investigation of energy efficiency of a commercial 5kW PEMFC fuel stack and on the possible directions to reduce the dimensions and weight of the PEMFC stack. Voltage dependencies were determined between voltage (U) -current (I) and current (I) – Power (P), as were the temperature distribution during the operation of the fuel cell stack and the peak current in humidifying the PEMFC stack by means of the SCU system. The energy demand of the 5kW fuel cell stack for the needs of the cooling system operation was determined (the so-called internal requirements of the fuel cell stack)

Streszczenie W pracy scharakteryzowano możliwości zastosowania niskotemperaturowych wodorowo-tlenowych ogniw paliwowych z protonowymienną membraną PEMFC w transporcie i lotnictwie. W pracy nacisk położono na zbadanie efektywności energetycznej komercyjnego stosu ogniw paliwowych PEMFC o mocy 5kW, a także możliwe kierunki zmniejszania jego masy i gabarytów. Wyznaczono zależności napięcie (U)– prąd (I) oraz prąd(I)-moc (P), rozkład temperatur podczas pracy stosu ogniw paliwowych, a także wielkości prądu w „piku” podczas nawilżania stosu PEMFC za pomocą układu SCU. Wyznaczono zapotrzebowanie energetyczne stosu ogniw paliwowych 5kW na potrzeby pracy układu chłodzenia (tzw. potrzeby własne stosu ogniw paliwowych) (Określenie parametrów elektrycznych oraz efektywności chłodzenia powietrznego niskotemperaturowego stosu ogniw paliwowych PEMFC o mocy 5kW)

Słowa kluczowe: ogniwa paliwowe z protonowymienną membraną PEMFC, chłodzenie powietrzne, wodór, efektywność energetyczna
Keywords: proton exchange membrane fuel cell PEMFC, air cooling, hydrogen, energy efficiency

Introduction

Among the five advanced types of hydrogen–oxygen fuel cells, the Polymer Membrane Fuel Cells (PEMFC) are more and more used in portable, stationary, and transport energy systems in different economic sectors. Nowadays, in Poland and worldwide, more and more attention is paid to the application of fuel cells (FC) in construction of power units for electric vehicles of a different type, unmanned aerial vehicles, aircrafts, means of the maritime and land transport, including submarines, inspection facilities for seabed, etc. [1-3]. According to the Fuel Cell Industry Review, in 2016, for the first time the total power of the fuel cells supplied for transport sector (280MW) was higher than the power of energy generators with fuel cells designated for stationary solutions (200MW). The main reason was the introduction of hydrogen-oxygen PEMFCs for construction of power units fueling the cars of Toyota Mirai in Japan, California and to a lesser extent in Europe [4]. The range of cars amounts to approx. 500-700 km, while the charge procedure for hydrogen containers takes approx. 5-6 minutes. PEMFCs are also increasingly used for construction of power units powering trucks, buses, and forklift trucks [5]. Generators with fuel cells (FC) with lower power, ranging from 100W to approx. 5 kW, are also used to power electric engines in scooters, bicycles, wheelchairs, unmanned ground robots [6,7]. An outstanding example is the “Burgman” scooter – the first motorcycle using hydrogen-oxygen fuel cells for power system. The motorcycle was designed by Suzuki, while PEMFC system is provided by Intelligent Energy. The motorcycle received the European type-approval as the first fuel cell vehicle in the world. Another interesting product is the hydrogenoxygen PEMFC bicycle developed by the consortium composed of the following companies: i) Cycleurope, manufacturer of electric bikes, ii) Ventec, manufacturer of energy management systems in electrochemical energy sources, iii) Pragma, manufacturer of PEMFC. Another innovative solution of this venture was the application of chemical source of hydrogen in the form of hydrogen-rich solid chemical compound. Hydrogen to power is produced on-line in reactor when riding a bike, and directly consumed by PEMFC [8,9].

Currently, in aviation it is aimed at reducing: consumption of non-renewable fuels, CO₂ emissions and pollution, noise-reduction by using electric engines in power units and by using electrochemical energy sources to power the said engines. Hypothetically, three types of electrochemical sources (reservoirs) of the source of electrical power may be used in aviation electric drives: electrochemical batteries (EB), supercapacitors (SC) and fuel cells operated on hydrogen or alternative fuels. The most promising solution for construction of power systems is the application of fuel cells, which ensure the highest amount of energy generated by hydrogen-powered fuel cell stacks, stored under pressure in the composite, ultralight hydrogen container. With respect to fuel cells, electrochemical batteries, and even more supercapacitors, are characterized by unfavorable amounts of energy (amount of energy stored in the system with unit weight). However, their advantage is high power density (power generated/device weight). For this reason, it is necessary to consider the application of hybrid system with main fuel system in the form of fuel cells and energy reservoirs as EB or SC. It would support the main fuel system for starting, maneuvering or emergency purposes [10,11].

The application of PEMFCs in the technology of unmanned aerial vehicles (drones with airframe or rotary structures) causes the extension of flight time as compared to electrochemical batteries. The examples of models of UAV with hydrogen-oxygen fuel cells include: Helios (manufacturer of NASA), Antares DLR-H2 (Lange Aviation), Ion Tiger (U.S. Navy), Puma (AFRL, PTX, MCEL), Pterosoar (Horizon), HyFish (DLR, Germany), Mirador (DGA, France), Spider-Lion (AeroVironment, USA) or models by Georgia Inst of Technology (2006), KAIST (2007). The flight time for UAVs fuelled by gaseous hydrogen ranges from 0.2h to approx. 6h. It should be stressed that usually electrochemical energy sources used to power electric engines of UAV are built in the form of hybrid units: fuel cells cooperating with the battery [12,13]. PEMFCs have been successfully used to build power units with higher electric power. Powered sailplane Antares DRLH2, or Enfica-FC serves as good examples of the application of PEMFC in this industry [14,15]. In 2016 the German engineers from DRL carried out successful flight test of four-seater passenger aircraft HY4. During the 10-minute flight there were 2 pilots and 2 manikins on board the aircraft HY4. It uses hydrogen fuel cells enabling to achieve a cruising speed of 145 km/h and distance of 1500 km. The energy from the batteries is used during take-off and when landing [16].

Another important application area of fuel cells in aviation is Auxiliary Power Unit (APU). The fuel cells as auxiliary power units have already been successfully applied in large passenger aircrafts, such as Being, Airbus [17,18]. The advantage of APUs with fuel cells is the fact that they can operate independently, without main electric engines, and supply electricity during flight and ground operations for on-board infrastructure [19]. The FCs are also more and more used to fuel aviation ground support equipment (in English: Aviation ground support equipment – GSE). The article [20] presents the possibility of application of the PEMFC fuel cells stack with power of 5 kW to fuel mobile light tower at San Francisco International Airport, US. Other possibilities include the application of fuel cells for fueling buses transporting passengers or other aviation auxiliary vehicles [21].

However, despite the huge interest in fuel cells, the reference literature lacks systematic research regarding the energy efficiency of the PEMFC stacks. Target power generators, based on the PEMFC stacks, are equipped with many additional auxiliaries, which improve and support their proper operation. They include systems monitoring the humidification and dosage of gaseous reagents, cooling and temperature control system, hydrogen leakage detection sensors, start-up accessories and others.

Power generators with the PEMFC fuel cells in the scope of power ranging from 20 to 200 kW, dedicated for applications in land and air transport, are constructed on a modular basis. In order to construct the target generator, several units (fuel cell stacks) with lower electric power, electrically connected in serial or in parallel, are used to achieve the desired electrical parameters (such as voltage, current and electric power) tailored to the requirements of the fuelled electricity receivers. In the case of power modules from the FC with electric power of approx. 20- 50kW, these are usually the FC stacks with power of 5-10 kW [22]

The aim of this article is to determine characteristics of electrical parameters and energy efficiency of commercial PEMFC fuel cell stack 5kW as reference unit for construction of future power generators.

Experimental part

The subject of the research covers the gaseous hydrogen-fuelled low-temperature fuel cell stack with power of 5kW – H5000 (by Horizon, Singapore). The H5000 FC stack operates in the so-called “open-cathode” system (in English: open-cathode PEMFC stack). In this structure the elements of the PEMFC stack (single components of MEA and graphite bipolar plates) are arranged so that cathode areas to fuel with oxygen mainly from cooling air flow generated by the fan unit are available outside. The technical parameters of the PEMFC fuel cell stack selected for the research are included in the manufacturer’s specification [23]. Fig. 1 shows the photo of the PEMFC fuel cell stack H5000 (Horizon, Singapore). Fig. 1a shows the photo of the H5000 from the side of cooling fans. While Fig. 1b shows the photo of the H5000 PEMFC stack construction from the other side, i.e. air inlet to the cathode area.

Fig.1. Photo of an H5000 PEMFC stack
a) View from the side where 4 cooling fans were installed in the device
b) View from the side of the air inlet to the cathode area of the PEMFC stack

Fig.2. Diagram of the stationary set-up for investigation of the electrical parameters of the PEMFC stack H5000

Apparatus and Method of Measurements

Fig. 2 shows the diagram of stationary setup for determining current (I) and voltage (U) characteristics of the tested H5000 PEMFC stack. The H5000 PEMFC stack was fuelled with hydrogen cylinder (Air Liquid ALPHAGAZ 1 H2 – purity 5.0). The dry gaseous hydrogen under pressure pH₂=0.5bar +/-0.05 bar was supplied to the anode area of the H5000 PEMFC stack. The shut-off electronic valve NC (Normal Close) is mounted in the hydrogen supply valve. The operation of the electronic valve is regulated by the automatic control system of the operation of the FC stack. Directly at the outlet of unreacted hydrogen from the FC stack, the electronic valve NC “purge” is mounted, whose operation is controlled by the PEMFC stack controller. The “purge” electronic valve is used to periodically rinse the anode area of the fuel cell stack. The cooling system is a very important element of the tested PEMFC stack. In this solution it is a system of 4 fans, with power of 100 W each. The PEMFC stack operation controller is an integral component of the tested PEMFC stack. The PEMFC stack is humidified using the SCU (Short Circuit Unit) causing short periodical circuits of the PEMFC stack.

The H5000 PEMFC stack controller was powered from the external direct current power supply at voltage of 24VDC.

Methodology of Electrical and Thermal Measurements of the H5000 FC Stack

The characteristics of current (I) – voltage (U) and current (I) – power (P) were determined using the electronic load Chroma 63202 (2600W/0-50A/0-600V) enabling to load the stack up to the power of 2.5 kW. In order to increase the load of the tested H5000 PEMFC stack up to the volume exceeding nominal electric power of 5 kW, the electrical system was retrofitted with resistive loads. The volume of current (I) and voltage (U) from the PEMFC stack was measured using the multimeters. The value of current was measured on the shunt 150A/60mV using the Agilent 34411A multimeter, while the voltage U was measured using the Agilent 34410A multimeter. The use of transducers LEM 15 enabled to determine the current and electric power collected by the stack control system, predominantly for the purpose of cooling the H5000 stack. The data was archived by the measuring card LabJack U3- HV on the PC. During the electrical tests of the H5000 PEMFC stack with permanent or variable load, the temperature distribution along the entire length of the fuel cell stack was measured. The measurement of the FC stack temperature distribution field from the side of air inlet was carried out using the infrared thermal camera by NEC Thermo Tracer H2640 equipped with the bolometer infrared detector, resolution of 640×480 pixels, enabling the registration of temperature with a resolution of 0.03°C, and the emissivity value was at the level of e=0.95. The thermographs were analysed using the Thermography Studio Professional software.

Flow Tests and Characteristics of the H5000 FC Stack Cooling System

The analysis of air inflow speed distribution to inlet channels in bipolar plates in the H5000 FC stack was carried out using the combined, three-channel thermoanemometric sensor TURBULENCE METER type ATM 94 (Fig.3).

The used thermoanemometric sensor (Fig.3) enabling the measurement of the absolute speed and its components oriented in a classic, Cartesian coordinate system, is composed of three active elements in the form of tungsten fibres, 5μm thick and approx. 2mm long, spanning between the brackets. The brackets not only serve as mechanical attachment of fibre, but also electrical connection with the cable connector. Individual sensor fibres are placed so that they create the edge of the cube, whose axis overlaps with the sensor axis. The described sensor cooperated with the measuring card A/C PC LabCard PLC 814 installed on the PC, along with the measurement software developed by the Laboratory of Flow Metrology of the Strata Mechanics Research Institute of the Polish Academy of Sciences in Krakow. Before the use, the sensor was calibrated at the Calibration Laboratory for Ventilation Measuring Instruments of the Strata Mechanics Research Institute of the Polish Academy of Sciences in Krakow, holding the national accreditation confirmed by the certificate No. AP 118.

Fig.3. ATM 94 three-fibre sensor: enlarged view and construction schema

The flow characteristics of axial fans installed in the cooling system of the H5000 PEMFC stack were determined at the test stand prepared according to the standard PN-EN ISO 5801:2008. Test stand of the C type – measuring pipeline is installed on the suction side of the tested fan (6). The scheme of test stand is presented in Fig. 4

Fig.4. Set-up for measurement of the flow characteristics of the fans

According to the numbering in Fig. 4 the louvre damper responsible for flow throttling (1) is installed at the inlet, the next element is the thermo-hygrometer model HD48T by Delta Ohm (2). In order to harmonize the speed profile, the stream straightener (3) is installed behind the damper. The measuring set composed of the airflow element by STRA Dwyer (4), averaging the dynamic pressure profile in the cable, and differential pressure transmitter by Halstrup Walcher P26 was used for airflow measurements. At a distance of 3 averages before the inlet to the rotor, the pressure signal at suction was introduced to the pressure gauge VPT-100 by Voltcraft. The tested fans (6) are powered from the laboratory power supply (8). The rotations were measured using the optoelectronic tachometer (7).

Research Results Characteristics of Performance Parameters of the H5000 PEMFC stack

The main factor determining the suitability of the PEMFC stack as the basic component for construction of power generators is the possible electric power for supplying to the power unit or other electricity receivers. Each galvanic cell, and thus the fuel cell stack, is characterized by maximum power points, to which the pair of parameters correspond (voltage-current)_max. For certain types of galvanic cells, this point may not be visible on curve power (P) – current (I) (it often happens for primary cells), since its location is beyond the useful voltages or currents. In the case of fuel cell generators, the location of the point (voltage-current)_maxis usually clearly marked and plays an important role during the generator operation. When the fuel cell stack is being gradually loaded, it results in self-adjusting change of operating conditions (voltage drop, current increase, increase of hydrogen fuel consumption), which causes an appropriate increase of the power from the stack according to the growing demand. This mechanism will work effectively until it reaches the maximum power point, after exceeding this point, the power supplied by the cell stack will rapidly decrease, despite the increasing current load [24]. In the case of uneven demand for power in time, the often-used solution aiming at balancing the power is to use the assistance system in the form of peak batteries or supercapacitors [25, 26]

Fig. 5 shows the determined dependencies voltage (U) – current (I) and current (I) – power (P) for the tested H5000 FC stack in the scope of nominal power declared by the manufacturer..

Fig.5. Voltage (U) – current (I) and power (P) – current (I) dependencies determined for the H5000 PEMFC stack in the required operation conditions (nominal power)

Based on the presented in Fig. 5 dependencies between voltage (U) – current (I) and current (I) – power (P), it can be stated that for the tested PEMFC stack, the power of the PEMFC stack gradually increases up to the 5kW during the increase of current load. It is the value declared by the manufacturer as nominal power.

On the basis of this dependence, it may be concluded that the tested H5000 PEMFC stack has not achieved the power point Pmax yet, where after exceeding, there is usually the sharp decrease of electric power. In the initial part of the characteristics, voltage (U) – current (I) for the current ranging from 1 to approx. 10A, a significant voltage drop from the highest value of U_ocv = 114V (OCV Open Circuit Voltage) to the value of U = 85V can be observed. This fact is related to activation losses dominant within this scope of the FC operation. Within the research scope of the applied electric loads, the power of the H5000 device does not exceed 1kW, and the operating temperature approximates the ambient temperature (20^oC).

Further increase of electrical load of the H5000 PEMFC stack with current of 10A-70A causes a linear voltage drop resulting from ohm losses dominant within this part of characteristics of U=f(I). The nominal power point of 5kW declared by the manufacturer of the FC stack corresponds to the pair of U = 72V, I = 70A. The received results are compliant with the nominal electrical parameters of the H5000 PEMFC stack declared by the manufacturer [23].

However, for stationary, transport or aviation applications of the fuel cells for construction of power generators, the important factor is the information whether or not it is possible to increase power in case of intermittent demand by the power unit [27,28].

Fig.6 shows another dependence between voltage (U)- current (I) and power (P)-current (I) under expanded range of electrical load, i.e. above the nominal power of 5kW.

Fig.6. Voltage (U) – current (I) and power (P) – current (I) dependencies determined for the H5000 PEMFC stack under expanded range of electrical load

As it follows from the determined dependencies (Fig. 6) between voltage (U) – current (I) and power (P) – current (I), the H5000 fuel cell stack achieved the power of approx. 6 kW. Similarly to the previous dependencies (Fig. 5), the tested generator also failed to achieve the maximum power point P_max. These results confirm that the declared and achieved power value of 5kW, by the tested power source, is not the highest peak value for this device. The tested H5000 FC stack has some power reserve. Due to the durability of the PEMFC and repetitive performance parameters of the tested stack, the team determined the critical power value not higher than 6kW.

Apart from electricity, the waste heat and water are the output of the PEMFC stack. In practice, the efficiency of the PEMFC stacks amounts to approx. 50%. Thus, in the case of the tested device, in addition to electricity of 5 kW, also the heat of 5 kW was received, which must be removed from fuel cells in this H5000 PEMFC stack. The most commonly used methods of waste heat collection include: a) air cooling using fans, b) cooling using the liquid medium, c) cooling using the passive elements [29-31].

The selection of cooling technology depends on the volume of the power generated by the PEMFC stack, its application (portable, stationary) as well as structure of bipolar plates combining the PEMFCs in the stack. The problem of cooling the PEMFC stack is a complex issue having a very large impact on the energy efficiency of the device. Too low efficiency of the cooling system may cause local and rapid increase of temperature along the entire length of the FC stack, and thus, may cause rapid loss of humidity by nafion polymer membrane, and in extreme cases, may result in thermal damage to a single MEA (in English: Membrane Electrode Assembly; MEA, name refers to a single PEM fuel cell). The occurrence of temperature gradient in the FC stack may also generate thermomechanical stresses in bipolar plates causing their gradual degradation [32].

On the other hand, too high cooling air flow rate causes: a) excessive removal of humidity from cathode area, which may result in drying of MEA and increase of internal resistance of the FC, and b) supercooling, which prevents the achievement of optimal operating temperature by the FC stack, decreasing the electrochemical reaction time, resulting in the decrease of efficiency of the whole FC. Thus, the cooling system is an integral element of power generator with the PEMFC. This fact makes it necessary to determine the energy consumption by the FC stack for powering and controlling the cooling system, in other words, for the so-called FC stack own purposes. This factor determines the need to select flow machines (fans) operating under optimized conditions, i.e. with high efficiency. The authors determined the demand for electricity for powering axial cooling fans installed on commercial H5000 FC stack. It should be stressed that cooling fans are responsible for main consumption of electricity (the so-called own purposes) of the fuel cell stack during operation.

Fig. 7 shows the dependence of variation in electric power N_el necessary to produce air flow by the cooling fan. The tests were performed for various supply voltages ranging from 8 to 24 V.

Fig.7. Dependence of variation in electric power N_el on the air stream produced by the fan for various supply voltages (range: 8‒24 V).

As it follows from Fig. 7, the highest electric power consumption N_el amounting to approx. 100 W, at a voltage of 24V was registered for the air flow ranging from 0.02 to 0.07 m³/s, while the increase of the cooling air flow produced by the fan causes insignificant decrease of the consumed power N_el to the level of approx. 70-80 W. Similar dependencies can be observed for lower supply voltages. The decreased demand for electric power Nel for powering the fan, along with the increase of air flow, is caused by lower increase of total pressure ∆p and an increase of efficiency, which is presented in Fig. 8 and 9. Fig. 8 illustrates the dependence of total pressure increase ∆p on air flow. The tests were performed for supply voltages ranging from 8 to 24 V

Fig. 8 presents the characteristics of increase of ∆p in efficiency function for different supply voltages of the fan motor. These characteristics compared to the characteristics of pressure losses in the FC stack are used to determine the fan operation point. Optimally, these points should fall within the areas of the highest efficiency of fans. Fig. 9 shows the determined dependence of electrical efficiency of a single axial fan in the function of air flow. By analysing the diagrams from Fig. 8 and Fig. 9 as well as flow characteristics of the PEMFC stack (pressure losses), the air distribution control algorithm, which is optimal in terms of energy, can be determined, which leads to the minimization of the cell own purposes

Fig.8. Dependence of total pressure ∆p on air flow at different supply voltages for a single fan.

Fig.9. Dependence of the total performance of the electric fan on its efficiency

Based on the presented flow characteristics of the fan, it can be stated that the cooling system composed of 4 fans, is characterized by moderate efficiency. For optimized operating conditions for a single device, it amounts to approx. 40-45% depending on the supply conditions. The air flow produced by axial fans through cathode channels in bipolar plates is mainly used to cool the FC stack, but it is also used to fuel the cathode of individual cells in the FC stack with oxygen. There are 120 single MEA units in the H5000 FC stack. Fig. 10 shows the dependence of the air inflow speed to the FC stack on the electric power consumed by the 4 axial fans cooling FC.

Fig.10. Dependence of the air inflow speed to the cathodic areas of the OP stack on changes in the amount of electric power consumed by the system of 4 fans.

As it follows from Fig. 10, the air inflow speed to the cathode areas of the FC stack, as expected, depends on the power consumed by the system of 4 fans and does not depend whether or not the air inflow speed to the PEMFC increases or decreases during the measurements. This fact is particularly important in dynamically changing conditions of the power generator operation with the FC, where it will be necessary to dose the same amount of air in order to maintain the performance parameter repetitiveness.

Another proper operation test of the cooling system was the analysis of homogeneity of air flow speed distribution to the cathode areas in the FC stack.

The flow qualitative test (comparison of the air inflow distribution) using the thermoanemometric probe was carried out within ¼ of the air inflow area to the H5000 FC stack. During the tests, the axial symmetry of flow phenomena in the tested H5000 device was assumed. Fig. 11 presents the photograph of the H5000 PEMFC stack with distribution of measuring series.

Fig.11. Analysis of homogeneity of air inflow speed to the cathodic areas of the PEMFC stack H5000.

In turn, Fig. 12 presents the determined changes in air inflow speed to the H5000 fuel cell stack according to the adopted measuring series. For the distance, the value “0” was assumed on the edge of the FC stack.

Fig.12. Dependence of changes in air inflow speed on the distance from the edge of the PEMFC stack H5000.

Based on the conducted measurements, it is stated that there is the differentiation of air inflow speed supplied by the fans of the PF stack. The air is supplied at different speed to individual single PEMFC in the FC stack. However, on the basis of electrical characteristics, a negative impact of this phenomenon is not stated in the case of free position of the device on stationary stand. If the device is installed in electric vehicle or aircraft airframe, this factor may have an impact limiting the FC stack to achieve the highest electrical parameters.

During the operation of the FC stack with variable electrical load, in addition to electricity, also heat is generated on the internal resistance of the FC as a result of electrical current flow. In order to keep the optimal operating temperature below 65^oC in the FC stack, it is cooled by means of the forced (by 4 fans) cooling air flow through cathode channels in bipolar plates. The temperature change distribution on the surface of the FC stack from the side of air inflow to cathode channels during the operation was recorded using the thermographic camera Fig.13 a-b.

Fig. 13b shows the temperature distribution profile on line 1 (Fig. 13 a) over the entire width of the FC stack during the maximum electrical load of the stack during the measurements. The even temperature distribution with slight increase of temperature of 2-4K at the ends of the FC stack can be observed both on the thermograph (Fig. 13a) and on line 1 (Fig. 13 b). It probably follows from minor unevenness of cooling air flow through cathode channels caused by 4 axial cooling fans positioned on the opposite side of the FC stack.

Fig.13a Image of temperature distribution recorded with a thermographic camera during H5000 PEMFC stack operation under load with a nominal power of 5kW Fig. 13b Temperature distribution profile on line 1 (from Fig. 13a)

An important issue also pertains to the analysis of humidification level of each 120 single PEMFC membrane in the H5000 FC stack. During the operation of the FC stack with current load, a crucial issue is the control of polymer membrane humidification level, so that its internal resistance is not increased due to the excessive drying. The H5000 FC stack can be humidified by retrofitting the reagent dosing system in membrane humidification devices. This solution is the most frequently applied for the FC stacks constructed in the so-called closed cathode system. On the other hand, in the case of the tested structure of the FC stack, with the so-called open cathode, a special SCU (Short Circuit Unit) is used for self-humidification of polymer membranes of the FC stack. As a result of short circuits of the FC due to electrochemical reaction, there is the excess water resulted from the increased flow of short-circuit current, which humidifies polymer membranes. A positive impact of the SCU on the H5000 stack performance parameters is presented in Fig. 14, where the compared characteristics voltage (U) – current (I) and power (P) – current (I) with the SCU switched on or off are presented.

Based on the recorded characteristics of voltage (U) – current (I) and power (P) – current (I) for the FC stack, with the SCU switched on or off respectively, it may be concluded that the fact that the SCU was switched on, it caused the improvement of both value of voltage and power generated by the FC stack. This fact is directly connected with the drop of its internal resistance as a result of humidification of the PEMFC membrane.

Fig.14. Comparison of the electrical characteristics U = f (I), P = f (I) of the PEMFC operating with the SCU (Short Circuit Unit) switched on or off.

One of the factors guaranteeing the continuity of humidification of the PEMFC membrane in the fuel cell stack by means of the SCU is the control of the SCU switching frequency. For the functioning of the control system of the FC stack performance parameters and other components of the power system, it is important to know the value of short-circuit current during switching on the SCU and the FC stack voltage reaction during and after the short circuit.

The knowledge of the above, fast-changing and short-term changes of current and voltage during the SCU short circuits should facilitate the selection of protective and receiving devices cooperating with the fuel cell stack, such as voltage converter, or necessity to retrofit the power system in supercapacitor.

For this purpose, during the short circuit caused by the SCU, the voltage drop at the output of the H5000 FC stack (comprising 120 fuel cells) as well as short-circuit current on the shunt 1000A/100mV were measured using the 2- channel oscilloscope RIGOL DS1062CA.

Fig. 15 shows a single SCU short circuit at load current of 36A, the duration of which amounts to approx. 55ms (curve 1 – short-circuit current measured on the shunt) and the corresponding cell voltage reaction (curve 2 – voltage on the fuel cell stack, multiplier x10).

Fig.15. Single SCU short circuit at load current 36A; duration of short circuit: 55 ms; curve 1: short-circuit current, curve 2: voltage of stack (multiplier x10).

Based on the performed weight analysis of individual components of the H5000 device, it is concluded that total weight of fans amounts to 3.3 kg, and the FC controller weight with electric cables amounts to approx. 3 kg. It should be stressed that the exchange of fans into the unit with higher efficiency and less weight when designing the device for portable solutions is not the serious problem for the present solutions in the technology. Another action aiming at reducing the FC stack weight is the exchange of casing made of aluminium (weight of approx. 6 kg) with the casing made of composite coat material or other plastic. The greatest technological problem is to replace bipolar graphite plates (estimated weight of 121 bipolar plates amounts to approx. 9kg, and the thickness of a single plate h_c amounts to 4mm, i.e. 121 plates in the H5000 FC stack amounts to 48 cm, components MEA 12 cm, each of 120 has 1mm), and the length of the whole H5000 FC stack with the casing amounts to approx. 65 cm. On the basis of the analysis of reference literature, it can be stated that the application of metal bipolar plates could enable reduction of the weight to approx. 4 kg, and reduction of the length of the whole FC stack by at least 30%.

Fig.16. A single bipolar plate made of graphite.

Summary

This article presents the research results of energy efficiency of the H5000 PEM fuel cell stack, which is air-cooled. The basic electrical parameters of the H5000 FC stack are determined in order to not only define the nominal power of the H5000 stack, but also to define the possibility to increase such power as a result of sharp increase in demand for electric power of the receiver. The characteristic parameters of the cooling system composed of 4 cooling fans are also determined. As it follows from the presented research, the energy efficiency of the applied fans ranges from 40 to 45%. A moderate efficiency of the FC cooling unit may be the reason for the increase of the fuel cell generator’s own purposes. The next analysed parameter was the speed distribution of air inflow to cathode openings in bipolar plates of the FC stack. Based on the presented research, the heterogeneous air inflow is determined, which does not materially affect the electrical parameters of the FC stack operating on a stationary stand in the ambient air cooling conditions. The PEMFC stacks for structures, in which metal bipolar plates enabling the reduction of weight and size are applied, will be the future solution for portable applications.

Acknowledgments The research presented in this paper was financed as a project PBS3/A6/24/2015 AOS-H2 the Applied Research Programme (PBS) of the National Centre for Research and Development (NCBIR), Poland, in the years 2015–18. Some of the measurements were performed using the research infrastructure of the AGH Centre of Energy.

REFERENCES

[1] Garche J, Encylopedia of Electrochemical Power Sources Elsevier 2009
[2] Wu H, A review of recent development: Transport and performance modelling of PEM fuel cells, Applied Energy, 1 (2016) 81-106
[3] Broussely M, Pistoia G, Industrial Application of Batteries. From Cars to Aerospace and Energy Storage, Elsevier, Amsterdam, 2007
[4] The Fuel Cell Industry Review 2016 http://www.fuelcellindustryreview.com (15.07.2017)
[5] Miranda P.E, Carreira E.S., Icardi U.A., Nunes G.S.Brazilian hybrid electric-hydrogen fuel cell bus: Improved on-board energy management system, International Journal of Hydrogen Energy, 42 (2017) 13949-13959
[6] Mellino S, Petrillo A, Cigolotti V, Autorino C, Ulgiati S, A Life Cycle Assessment of lithium battery and hydrogen-FC powered electric bicycles: Searching for cleaner solutions to urban mobility, International Journal of Hydrogen Energy, 42, (2017), 1830-1840
[7] Rajashekara K, Hybrid and fuel cell system for transportation, Electrical Review, 80 (2004) 1033-1039
[8] London police trialing Suzuki scooters with Intelligent Energy units, Fuel Cell Bulletin, 9 (2017) 3-4
[9] French group launch fuel cell powered bike to boost emobility, Fuel Cell Bulletin, 7 (2013) 3-4
[10] Pratt J., Klebanoff L., Ramos K., Akhil A. A., Curgus D.B., Shenkman B.L., Proton exchange membrane fuel cells for electrical power generation on-board commercial airplanes, Applied Energy 101 (2013) 776-796
[11] Bicer Y, Dincer I, Life cycle evaluation of hydrogen and other potential fuels for aircrafts, International Journal of Hydrogen Energy, 42 (2017) 10722-10738
[12] Donateo T., Ficarella A., Spedicato L., Arista A., Ferraro M. A new approach to calculating endurance in electric flight and comparing fuel cells and batteries, Applied Energy 187 (2017) 807–819
[13] MMC next-gen fuel cell powered drone with telecom upgrade Fuel Cells Bulletin, 3 (2017) 5
[14] Intelligent Energy fuel cells significantly boost drone flight time, Fuel Cells Bulletin, 1 (2016) 4
[15] Fuel cells power first take off for DLR’s Antares aircraft, Fuel Cell Bulletin, 9 (2009) 3
[16] Romeo G., Borello F., Correa G., Cestino E., ENFICA-FC: Design of transport aircraft powered by fuel cell & flight test of zero emission 2-seater aircraft powered by fuel cells fueled by hydrogen, International Journal of Hydrogen Energy, 38 (2013) 469-479
[17] HySA Systems, SA National Aerospace Centre, Airbus study fuel cells for commercial airliners, Fuel Cell Bulletin, 9 (2014) 1
[18] Boeing fuel cell plane in manned aviation first, Fuel Cell Bulletin, 4 (2008) 1
[19] Sharaf O.Z., Orhan M.F., An overwiew of fuel cell technology: Fundamentals and application, Rewenable and Sustainable Energy Reviews, 32 (2014) 810-853
[20] Multiquip launches lighting tower powered by Altergy fuel cell, Fuel Cells Bulletin, 2 (2011) 5
[21] Plug Power, FedEx project rolls out fuel cell airport tractors, Fuel Cells Bulletin, 5 (2015) 2-3
[22] Fontela A., Soria A., Mielgo J., Sierra J.F., Blas J., Gauchina L., Martinez M., Airport electric vehicle powered by fuel cell, Journal of Power Sources, 10 (2007) 184-193
[23] https://www.horizonfuelcell.com/h-series-stacks (11.10.2017)
[24] Dudek M., Celowski P., Lis B., Raźniak A., Dudek P., Laboratoryjny generator energii elektrycznej o mocy 360W zawierający niskotemperaturowy stos ogniw paliwowych PEMFC chłodzony za pomocą medium ciekłego, Electrical Review, 92 (2016) 235–242
[25] Lü X., Miao X., Liu W., Lü J., Extension control strategy of a single converter for hybrid PEMFC/battery power source, Applied Thermal Engineering, 128 (2018) 887-897
[26] Lopez G.L., Rodriguez R.S., Alvarado V.M., Gomez-Aguilar J.F., Mota J.E., Snadoval C., Hybrid PEMFC-supercapacitor system: Modeling and energy management in energetic macroscopic representation, Applied Energy, 205 (2017) 1478-1494
[27] Sami B.S, Sihem N., Zafar B., Adnane Ch., Performance study and efficiency improvement of Hybrid Electric System dedicated to transport application, International Journal of Hydrogen Energy, 27 (2017) 12777-12789
[28] Saadi A., Becherif M., Hissel D., Ramadan H.S., Dynamic modeling and experimental analysis of PEMFCs: A comparative study, International Hydrogen Energy, 42 (2017) 1544-1557
[29] Chen, Y.Ch., Huang K.P., Yan W.M., Lai M.P., Yang C.C., Development and performance diagnosis of a high power aircooled PEMFC stack, International Journal of Hydrogen Energy, 41 (2016) 11784-11793
[30] Fly A., Thirng R., A comparison of evaporative and liquid cooling methods for fuel cell vehicles, International Journal of Hydrogen Energy, 41 (2016) 14217-14229
[31] Özden E., Tolj I., Barbir F., Designing heat exchanger with spatially variable surface area for passive cooling of PEM fuel cell, Applied Thermal Engineering, 51(2013) 1339-1344
[32] Zhang G., Kandlikar S.G., A critical review of cooling techniques in proton exchange membrane fuel cell stacks, International Journal of Hydrogen Energy, 37 (2012) 2412-2429

Authors: dr hab. inż. Magdalena Dudek AGH Akademia Górniczo – Hutnicza im. Stanisława Staszica w Krakowie, Wydział Energetyki i Paliw , al. A. Mickiewicza 30, 30-059 Kraków, E-mail: potoczek@agh.edu,pl; dr inż. Piotr Dudek, AGH-Akademia Górniczo – Hutnicza im. Stanisława Staszica w Krakowie, Wydział Inżynierii Mechanicznej i Robotyki, al. A. Mickiewicza 30, 30-059 Kraków, E-mial: pdudek@agh.edu.pl; mgr inż. Wojciech Kalawa, AGH Akademia Górniczo – Hutnicza im. Stanisława Staszica w Krakowie, Wydział Energetyki i Paliw, al. A. Mickiewicza 30, 30-059 Kraków e-mail: kalawa@agh.edu.pl; dr inż. Andrzej Raźniak AGH- Akademia Górniczo – Hutnicza im. Stanisława Staszica w Krakowie, Wydział Energetyki i Paliw, al. A. Mickiewicza 30, 30-059 Kraków, E-mial: razniak@agh.edu.pl; dr inż. Tomasz Siwek, AGHAkademia Górniczo-Hutnicza im. Stanisława Staszica w Krakowie, Wydział Energetyki i Paliw, al. Mickiewicza 30, 30-059 Kraków, Emial: siwek@agh.edu.pl;

Source & Publisher Item Identifier: PRZEGLĄD ELEKTROTECHNICZNY, ISSN 0033-2097, R. 94 NR 4/2018. doi:10.15199/48.2018.04.34