Next Article in Journal
Analyzing Potential Failures and Effects in a Pilot-Scale Biomass Preprocessing Facility for Improved Reliability
Next Article in Special Issue
Novel Hierarchical Energy Management System for Enhanced Black Start Capabilities at Distribution and Transmission Networks
Previous Article in Journal
A Long-Term Power Supply Risk Evaluation Method for China Regional Power System Based on Probabilistic Production Simulation
Previous Article in Special Issue
Reed Switch Overcurrent Protection: New Approach to Design
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Analysis of Power System Electromagnetic Transients Using the Finite Element Technique

Department of Electrical Power Engineering, Faculty of Electrical Engineering, Mechanical Engineering and Naval Architecture, University of Split, R. Boskovica 32, 21000 Split, Croatia
*
Author to whom correspondence should be addressed.
Energies 2024, 17(11), 2517; https://doi.org/10.3390/en17112517
Submission received: 25 March 2024 / Revised: 10 May 2024 / Accepted: 22 May 2024 / Published: 23 May 2024
(This article belongs to the Special Issue Energy, Electrical and Power Engineering 2024)

Abstract

:
In this paper, a numerical model for the analysis of electromagnetic transients in a power system, based on the finite element technique, was developed. The simplicity of the finite element technique is manifested in the fact that the problem of solving any mathematically described phenomena in an area is reduced to solving those same phenomena in a small part of that area, i.e., finite elements. Based on appropriate mathematical models, numerical models of the synchronous generator and other parts of the power system have been developed. System of differential equations of each power system element have been included in the numerical model in such a way that we employ numerical integration of the differential equations using the generalized trapezoidal rule ( ϑ -method) and reduce it to a system of algebraic equations in each time step. In a practical sense, this method enables a very simple solution to the problem of a real power system, i.e., a system with many generators, transformers, transmission lines and other elements. The accuracy of the developed numerical model for the analysis of electromagnetic transients in a power system has been confirmed by comparing the calculated results with the results obtained using the EMTP software (version 3.3) package.

1. Introduction

Analysis of electromagnetic transients in power systems is essential for understanding how the power system behaves during sudden changes. Several software options are available for this purpose [1,2,3]. Commonly used programs like EMTP [4,5], EMTDC [6], SPICE [7], and MATLAB–SIMULINK [8,9] employ electromagnetic time domain transient simulation techniques, utilizing methods such as mesh, state variables, and modified nodal analysis (MNA). Among these, the Electromagnetic Transients Program (EMTP) is widely preferred in scientific and engineering communities since it tackles a broad array of electromagnetic transient problems using a blend of the trapezoidal rule, the method of characteristics, and MNA, thus enabling accurate and robust simulation. This program package will thus serve as an excellent generator of benchmark results which will be used to compare results obtained by the numerical model proposed in this paper.
The proposed numerical model for the analysis of electromagnetic transients in the power system is based on the finite element technique (FET). This technique has been widely used for the solution of a great number of engineering problems in structural analysis and various continuum field problems [10,11]. The simplicity of the FET is manifested in the fact that the problem of solving any mathematically described phenomena in an area is reduced to solving this phenomenon in a small part of that area, which is, in the FET terminology, referred to as a finite element. The FET procedure has already been successfully applied to power system fault studies [12] as well as used for numerically solving various transmission line problems [13,14] and power system transients stability problems [15]. For example, in [15], an advanced FET-based program for numerical analysis of power system transient stability has been developed where variables of each finite element are defined as dynamic phasors (frequency domain) during the entire time-domain simulation.
As for the proposed numerical model for the analysis of electromagnetic transients in the power system, the mathematical model of the synchronous generator in the time domain has been developed which is defined and treated as a finite element with three local nodes. This original model includes an already developed numerical model of the excitation system as well as the model of the turbine governing system. Furthermore, based on the appropriate mathematical models, an originally developed numerical models of other parts of the power system (transformers, transmission lines, equivalent network and load) have been developed. They are each treated as separate finite elements in the computation process. The system of differential equations which encompasses all these elements is in the model reduced to a system of algebraic equations using numerical integration, more specifically, the generalized trapezoidal rule ( ϑ -method). Although the integration parameter ϑ is generally defined in the range from 0 to 1, there is restriction of its value regarding the stability properties of the time-stepping scheme. It can be noted here, for example, that various values of ϑ yield different finite difference schemes. These include, for example, the forward difference (Euler) scheme ( ϑ = 0 ), central difference (Crank–Nicolson) method ( ϑ = 1 / 2 ) and Galerkin scheme ( ϑ = 2 / 3 ) and backward difference scheme ( ϑ = 1 ) to name a few. The influence of ϑ on the stability of one-step time integration methods for the first-order initial value problems is discussed in [11,16], where the stability behavior is investigated by examining the free response of the finite element model. It was found that methods for which 0.5 ϑ 1 are unconditionally stable. Specifically, in this paper and the proposed model, the Galerkin scheme ( ϑ = 2 / 3 ) was used for time integration to maintain high stability of the results. In summary, by using the generalized trapezoidal rule, the system of differential equations is practically reduced to a system of algebraic equations in each time step. Depending on the finite element type, each finite element comprises a certain number of local nodes. The number of local nodes of each finite element is equal to the number of the physical connections of the element to the rest of the power system. In a practical sense, this method enables a very simple solution to the problem of a real power system, i.e., a system with many generators, transformers, transmission lines and other elements.
It is also important to note that the simultaneous simulation of both electromagnetic and electromechanical transients has always posed a considerable challenge to scientists and engineers. Previous efforts to tackle this problem are reported in [17,18,19,20,21,22,23]. Unlike EMTP, the FET procedure encompasses an assembly procedure which based on the connectivity matrix practically connects local nodes of finite elements to global nodes of the observed electric power system. For example, in the case of multiconductor transmission line (MTL) fast transient problems [14], finite element method analysis of a network with MTL and lumped elements is performed in such a way that a part of the network with MTL is subjected to the spatial discretisation into one-dimensional finite elements, whereas the electrical lumped elements are treated as separate finite elements. The corresponding connectivity matrix connects these elements into one global system.

2. Overview of the Developed Numerical Model

In this paper, we will demonstrate a successful application of the FET to modeling a complex electric power system during transient states. By applying the FET, the electric power system is divided into smaller parts, all of which are treated as separate finite elements. Each individual finite element is described by a set of ordinary differential equations which are solved numerically using the ϑ -method for a sufficiently small time interval. In such a way, all ordinary differential equations are practically transformed into algebraic equations which will be used to define the local system of equations for each finite element. The number of local nodes per finite element will be equal to the number of physical connections of that element to the electric power system. So for example, a synchronous generator finite element will have three local nodes, while the finite element of a three-phase transformer of a connection group Dyn5 will have seven local nodes.
As the most significant part of the power system, in the first step, a numerical model of the synchronous generator in the time domain has been developed, which is defined and treated as a finite element with three local nodes. Armature winding on stator, excitation and damping winding on a rotor are taken into account and are mutually coupled. The time-domain model of the turbine governor (Appendix A) as well as excitation system (Appendix B) have been already developed in [15] and incorporated into time-domain synchronous generator finite element model suitable for electromagnetic transient analysis. Single-phase and three-phase transformers are modeled as finite elements with four local nodes in the case of a single-phase transformer and a different amount of local nodes depending on the connection group in the case of a three-phase transformer. Single-phase and three-phase star-connected loads are modeled as finite elements with two local nodes and four local nodes for the single-phase and three-phase case, respectively. The equivalent network is modeled as a finite element with three local nodes. The complete theory concerning the synchronous generator, transformers, loads and equivalent network is originally developed in this paper. The time-domain model of a multi-conductor transmission line has been already developed in [14]. A three-phase transmission line as a finite element is represented by six local nodes where the local matrix and vector suitable for the assembling procedure with other power system elements have been defined in [14]. In this paper, we modeled only the basic elements of the power system. It is important to emphasize that it is possible to obtain a numerical model of any element of the power system. It is only important to know its mathematical model to define it as a finite element with a certain number of local nodes using the same procedure. The new approach enables, if necessary, a more detailed modeling of a particular part of the power grid, so that depending on the type of transient phenomena, different mathematical and therefore numerical models of a particular part of the power system can be defined.
Before the start of the transient analysis of the observed electric power system ( t = 0 ), it is necessary to initialize real initial conditions in the system itself. The initialization procedure is accomplished in a way that the real and reactive powers of generators are prescribed, whereas all other currents and voltages in the system are set to zero. The real powers of generators are prescribed in per-unit values via the mechanical turbine moments, whereas the reactive powers of generators are prescribed indirectly via the excitation regulator, i.e., the previously computed initial excitation voltage. Based on the known initial conditions and based on known (and previously derived) local systems of equations for all finite elements in the observed electric power system, assembling of these local systems into a global system of the power system is performed via the connectivity matrix which practically connects local nodes of finite elements to global nodes of the observed electric power system. This connectivity matrix is individual to each observed electric power system. This is standard FET assembling procedure [10,11]. To carry out assembling procedure the local systems of each finite must be set in the general form of
i l o c + = A l o c · φ l o c + + B l o c
where i l o c + is the vector of local nodal currents at the end of the integration interval, A l o c is the local matrix of the finite element, B l o c is the local vector of the finite element, and φ l o c + is the vector of local nodal potentials at the end of the integration interval.
By assembling all finite elements (parts of the power system), the following global system is obtained:
A g l o b · φ g l o b + = B g l o b
where A g l o b and B g l o b are the global matrix and vector which are formed by assembling all local matrices and vectors using the connectivity matrix, φ g l o b + is the vector of global potentials at the end of the integration interval.
By solving the global system of equations, unknown global potentials are obtained. After solving the global system of equations, global potentials are inserted back into finite elements to obtain local potentials again via the connectivity matrix. On the basis of known potentials of finite elements on the start and the end of the time interval, and the known currents at the beginning of the time interval, all unknown currents at the end of the time interval are computed along with all other additional variables of each finite element. The procedure is repeated for the next time increment bearing in mind that the quantities computed at the end of the previous time interval are now quantities at the beginning of the next time interval.
Based on the presented theory, a MATLAB program for electromagnetic transient analysis of power system has been developed. The accuracy of the computer program, based on the observed electric power system in test example, has been compared with solutions obtained by a well-known and established EMTP simulation tool, where a very good correlation has been demonstrated.

3. Power System Parts Modeled as Finite Elements

3.1. Synchronous Generator Modeled as a Finite Element

We will observe a three-phase synchronous generator with an armature winding on the stator, and a field winding (excitation winding) and a damping winding on the rotor. The damping winding is further divided into an equivalent winding in the d-axis and an equivalent winding in the q-axis. The spatial distribution of all windings on the rotor and stator is depicted in Figure 1.
Note that the proposed mathematical model of a synchronous generator will employ a system of natural coordinates ( a , b , c ) depicted on Figure 1.
According to Figure 1, we can write equations for all flux linkages ( Ψ ) in matrix form taking into account mutual inductances between different windings as
Ψ a Ψ b Ψ c Ψ f Ψ D Ψ Q = L a ( γ ) L a b ( γ ) L a c ( γ ) L a f ( γ ) L a D ( γ ) L a Q ( γ ) L b a ( γ ) L b ( γ ) L b c ( γ ) L b f ( γ ) L b D ( γ ) L b Q ( γ ) L c a ( γ ) L c b ( γ ) L c ( γ ) L c f ( γ ) L c D ( γ ) L c Q ( γ ) L f a ( γ ) L f b ( γ ) L f c ( γ ) L f ( γ ) L f D ( γ ) 0 L D a ( γ ) L D b ( γ ) L D c ( γ ) L D f ( γ ) L D ( γ ) 0 L Q a ( γ ) L Q b ( γ ) L Q c ( γ ) 0 0 L Q ( γ ) · i a i b i c i f i D i Q
where
  • Ψ represents the flux linkage of individual generator winding;
  • L represents the inductance of individual generator winding and mutual inductances between different generator windings;
  • γ is the electrical angle of the rotor position;
  • i is the current flowing through an individual winding;
  • a , b , c are indexes associated with the armature winding on the stator (see Figure 1);
  • f is an index associated with the excitation winding on the rotor (see Figure 1);
  • D , Q are indexes associated with the damping winding on the rotor (see Figure 1).
In the second step of deriving the mathematical model of the generator, we employ the well-known voltage equations of the stator and rotor windings in the natural coordinates and write them in matrix form as follows:
u a u b u c u f 0 0 = r a 0 0 0 0 0 0 r b 0 0 0 0 0 0 r c 0 0 0 0 0 0 r f 0 0 0 0 0 0 r D 0 0 0 0 0 0 r Q · i a i b i c i f i D i Q + d d t Ψ a Ψ b Ψ c Ψ f Ψ D Ψ Q
where r represents the resistance of the individual winding. Note that the voltages of equivalent damping windings D and Q are zero since the winding is short-circuited. The previous matrix equation can be written more concisely as follows:
u = R · i + d d t Ψ
In the finite element terminology, the number of local nodes of the finite element model of the synchronous generator will equal the number of physical connections of that finite element to the electric power system. In this case, there will be three local nodes depicted in Figure 2. As seen from the same figure, to each local node, a corresponding current will be joined as well as a corresponding voltage.
In addition to Equation (5), three additional swing equations are needed to obtain the local system of equations of the finite element model of a synchronous generator. These are needed since the matrix of inductances (3) is directly dependent upon the rotor position, i.e., the rotor movement:
ω r = d γ d t
T m · d ω r d t = m m m e l
s = 1 ω r = d δ d t
where
  • s is the synchronous generator rotor slip in p.u.;
  • ω r is the rotor angular frequency in p.u.;
  • δ is the rotor angle in r a d ;
  • T m is the mechanical time constant in r a d ;
  • m m is the per-unit (p.u.) mechanical moment of the generator;
  • m e l is the p.u. electrical moment of the generator.
As mentioned before, numerical integration will be employed to obtain simple algebraic equations of the generator instead of differential equations given by (5)–(8). These differential equations are firstly written implicitly and prepared for numerical integration over a sufficiently small time increment starting from t and ending at t + :
t t + ( d d t Ψ + R · i u ) · d t = 0
t t + ( d γ d t ω r ) · d t = 0
t t + ( T m · d ω r d t m m + m e l ) · d t = 0
t t + ( d δ d t s ) · d t = 0
Time integration of Equations (9)–(12) by the ϑ -method [10,11] yields the following equations:
Ψ + Ψ + R · ϑ · Δ t · i + + R · 1 ϑ · Δ t · i ϑ · Δ t · u + 1 ϑ · Δ t · u = 0
γ + = γ + ω r · Δ t
ω r + = ϑ · Δ t T m · m m + + 1 ϑ · Δ t T m · m m ϑ · Δ t T m · m e l + 1 ϑ · Δ t T m · m e l + ω r
δ + = δ + s + · ϑ · Δ t + s · ( 1 ϑ ) · Δ t
Note that the “+” sign denotes electromagnetic quantities at the end of the time interval. In Equation (13), substituting the vector of flux linkages with a inductance matrix multiplied by the current vector (see Equation (3)) both at the beginning and at the end of the time interval yields the following equation:
( R · ϑ · Δ t + L + ) · i + + ( R · ( 1 ϑ ) · Δ t L ) · i ϑ · Δ t · u + ( 1 ϑ ) · Δ t · u = 0
In order to obtain a complete local system of equations of the synchronous generator finite element model, one must rewrite the previous equation in a simpler form:
i + = D · u + + ( 1 ϑ 1 ) · D · u + E · i
where the auxiliary matrices are
D = F 1 · ϑ · Δ t
E = F 1 · ( L R · ( 1 ϑ ) · Δ t )
F = L + + R · ϑ · Δ t
It is important to note here that the generator was modeled as a finite element with three local nodes which means that it must be assembled with the rest of the power system elements as such. For that purpose, the matrix equation given by (18) must be written so that the stator quantities are separated from the rotor quantities. This will be accomplished using submatrices:
i 1 + i 2 + = D 11 D 12 D 21 D 22 · u 1 + u 2 + + ( 1 ϑ 1 ) · D 11 D 12 D 21 D 22 · u 1 u 2 + E 11 E 12 E 21 E 22 · i 1 i 2
where the newly introduced vectors are parts of previously defined vectors:
u 1 = u a u b u c i 1 = i a i b i c u 2 = u f 0 0 i 2 = i f i D i Q
The system of Equation (22) can be separated into two matrix equations of which only one is needed to be assembled with the rest of the power system elements and practically represents the local system of equations of the synchronous generator:
i 1 + = D 11 · u 1 + + D 12 · u 2 + + ( 1 ϑ 1 ) · D 11 · u 1 + ( 1 ϑ 1 ) · D 12 · u 2 + E 11 · i 1 + E 12 · i 2
Note that the vector u 2 + which practically represents the field winding voltage at the end of the time interval is determined on the basis of the voltage field regulator which is explained in Appendix B, while the mechanical moment of the generator m m is determined on the basis of the turbine regulator which is explained in Appendix A. The output mechanical power from the turbine regulator ( P m ) is equal in unit values to the mechanical moment of the generator ( m m ) since the rotor angular frequency ( ω r ) is also defined in p.u. and has the value 1 after the initialization procedure, i.e., at the time t = 0 .

3.2. Network Equivalent Modeled as a Finite Element

If a part of a power system is observed within one larger power system, it is necessary to separate the observed part of the power system from the remainder of the power system. This is accomplished by modeling the remaining power system using a network equivalent. Generally, to accurately model any network equivalent, it is necessary to know the following:
  • U n —nominal voltage of the network;
  • S k 3 —Three-phase subtransient short circuit power;
  • S k 1 —Single-phase subtransient short circuit power.
Based on this input data, one can determine the direct, inverse and zero sequence components of the equivalent reactance of the network:
X m d = X m i = U n 2 S k 3 Ω
X m 0 = 3 · U n 2 S k 1 2 · U n 2 S k 3 Ω
where
  • X m d is the direct sequence network reactance;
  • X m i is the inverse sequence network reactance;
  • X m 0 is the zero sequence network reactance.
Using the now-known symmetrical component reactances, one can determine the reactances, i.e., inductances in the natural a, b, c system:
X v = X a = X b = X c = X m 0 + 2 · X m d 3
X m = X a b = X a c = X b c = X m 0 X m d 3
It follows that that the inductances can be calculated by
L v = X v 2 · π · f
L m = X m 2 · π · f
This can be written in matrix form as
L = L v L m L m L m L v L m L m L m L v
The network equivalent as a finite element is depicted in Figure 3. The equivalent was modeled as a serial combination of electromotive forces and inductances taking the mutual inductance into account. Based on Figure 3, it is necessary to define a local system of equations for the finite element which consists of three local nodes.
The theoretical basis for forming the local system of equations for the network equivalent is the voltage equation:
φ u = L · d i d t
where the potential, current and voltage vectors are given by
φ = φ 1 ( t ) φ 2 ( t ) φ 3 ( t ) u = u a ( t ) u b ( t ) u c ( t ) i = i 1 ( t ) i 2 ( t ) i 3 ( t )
Using the ϑ -method, numerical integration of (32) is performed which yields the following result:
φ + · ϑ · Δ t + φ · ( 1 ϑ ) · Δ t u + · ϑ · Δ t u · ( 1 ϑ ) · Δ t L · ( i + i ) = 0
Rearranging (34) in a way that the current vectors at the end of the time interval are placed on the left-hand side, whereas all other quantities are placed on the right-hand side of the equation, the following is obtained:
i + = L 1 · ϑ · Δ t · φ + + L 1 · ( 1 ϑ ) · Δ t · φ L 1 · ϑ · Δ t · u + L 1 · ( 1 ϑ ) · Δ t · u + i
The system of equations provided in (35) represents the complete local system of equations of the network equivalent when modeled as a finite element.

3.3. Transformers Modeled as Finite Elements

The complete local system of equations for a single-phase transformer (without taking into account the transversal branch) will be derived starting from the equivalent scheme with the reduced secondary quantities, as depicted in Figure 4.
In the equivalent scheme of the single-phase transformer depicted in Figure 4, the following electromagnetic quantities are used:
  • U p ( t ) and U s ( t ) —voltages on the primary and secondary windings;
  • U s ( t ) —voltage on the secondary winding reduced to the primary side;
  • i p ( t ) and i s ( t ) —current through the primary winding and the current through the secondary winding reduced to the primary side;
  • R 1 and R 2 —resistance of the primary winding and the resistance of the secondary winding reduced to the primary side;
  • L σ 1 and L σ 2 —leakage inductance of the primary winding and the leakage inductance of the secondary winding reduced to the primary side.
The transformation ratio of the single-phase transformer equals
a = N 1 N 2
with N 1 and N 2 being the number of windings on the primary and secondary side. Using this ratio, one can describe the secondary side electromagnetic quantities reduced to the primary side as
U s ( t ) = a · U s
R 2 = a 2 · R 2
L σ 2 = a 2 · L σ 2
The basis for forming the local system of equations for a single-phase transformer is the following differential equation, which represents the mathematical model of the single-phase transformer according to (Figure 4):
U p ( t ) U s ( t ) = i p ( t ) · ( R 1 + R 2 ) + ( L σ 1 + L σ 2 ) · d i p ( t ) d t
Instead of solving this differential equation, one can choose a sufficiently small time increment over which numerical integration can be performed on (40) while maintaining a high degree of accuracy. This yields a simple algebraic equation instead of the differential equation. Specifically, time integration of Equation (40) is performed from some starting time t to an ending time t + :
t t + ( U p ( t ) U s ( t ) i p ( t ) · ( R 1 + R 2 ) ( L σ 1 + L σ 2 ) · d i p ( t ) d t ) d t = 0
Introducing Equation (37) into (41) and applying the ϑ -method on (41), the following equation is obtained after a certain amount of rearrangement to simplify the expression
i p + = G · U p + a · G · U s + + G · ( 1 ϑ 1 ) · U p a · G · ( 1 ϑ 1 ) · U s + H · i p
where
G = ϑ · Δ t ( R 1 + R 2 ) · ϑ · Δ t + ( L σ 1 + L σ 2 )
H = ( L σ 1 + L σ 2 ) ( R 1 + R 2 ) · ( 1 ϑ ) · Δ t ( R 1 + R 2 ) · ϑ · Δ t + ( L σ 1 + L σ 2 )
Quantities marked with a superscript “+” are quantities at the end of the time interval, while those without the “+” sign are those at the beginning of the time interval.
According to the finite element technique, the single-phase transformer can be modeled as one finite element with four local nodes (Figure 5). The electrical scheme on Figure 5 contains nodal currents and potentials. Note that the secondary side quantities are reduced to the primary side.
According to Figure 5 and the finite element technique terminology, voltage on the primary winding and the voltage on the secondary winding can be described by the following difference of local node potentials (this also extends to the quantities at the end of the time interval):
U p ( t ) = φ 1 ( t ) φ 2 ( t )
U s ( t ) = φ 3 ( t ) φ 4 ( t )
Furthermore local nodal currents of the finite element model are connected straightforwardly to the currents through the primary and secondary winding as follows:
i 1 = i p
i 2 = i p
i 3 = a · i 3 = a · i p
i 4 = a · i 4 = a · i p
Taking all of this into account and introducing Equations (45) to (50) into Equation (42), the following complete local system of equations of the four-node finite element of the single-phase transformer is obtained:
i 1 i 2 i 3 i 4 + = G · 1 1 a a 1 1 a a a a a 2 a 2 a a a 2 a 2 · φ 1 φ 2 φ 3 φ 4 + + G · ( 1 ϑ 1 ) · 1 1 a a 1 1 a a a a a 2 a 2 a a a 2 a 2 · φ 1 φ 2 φ 3 φ 4 + H · i 1 i 2 i 3 i 4
Now, in order to obtain the complete local system of equations for a three-phase transformer, one can assemble three single-phase transformers, depending on the connection group. The procedure of forming a complete local system of equations for a three-phase transformer will be elaborated on a Dy5 connected transformer which is often used in practice.
This Dy5 connected transformer will be modeled as a finite element with seven local nodes as depicted in Figure 6.
The three-phase Dy5 transformer is separated into three single-phase transformers as depicted in Figure 7.
According to Equation (51) and Figure 6 and Figure 7, the complete local system for the first, second and third pair of windings, respectively, can be written as follows:
i 11 i 21 i 71 i 41 + = G · 1 1 a a 1 1 a a a a a 2 a 2 a a a 2 a 2 · φ 1 φ 2 φ 7 φ 4 + + G · ( 1 ϑ 1 ) · 1 1 a a 1 1 a a a a a 2 a 2 a a a 2 a 2 · φ 1 φ 2 φ 7 φ 4 + H · i 11 i 21 i 71 i 41
i 22 i 32 i 72 i 52 + = G · 1 1 a a 1 1 a a a a a 2 a 2 a a a 2 a 2 · φ 2 φ 3 φ 7 φ 5 + + G · ( 1 ϑ 1 ) · 1 1 a a 1 1 a a a a a 2 a 2 a a a 2 a 2 · φ 2 φ 3 φ 7 φ 5 + H · i 22 i 32 i 72 i 52
i 33 i 13 i 73 i 63 + = G · 1 1 a a 1 1 a a a a a 2 a 2 a a a 2 a 2 · φ 3 φ 1 φ 7 φ 6 + + G · ( 1 ϑ 1 ) · 1 1 a a 1 1 a a a a a 2 a 2 a a a 2 a 2 · φ 3 φ 1 φ 7 φ 6 + H · i 33 i 13 i 73 i 63
According to Figure 6 and Figure 7, the following relations between nodal currents of the three-phase transformer relative to currents of local nodes of each of the winding pairs can be defined as follows:
i 1 + = i 11 + + i 13 + i 2 + = i 21 + + i 22 + i 3 + = i 32 + + i 33 + i 4 + = i 41 + i 5 + = i 52 + i 6 + = i 63 + i 7 + = i 71 + + i 72 + + i 73 +
So, by adding the contributions of all winding pairs, one can obtain a complete local system of equations of the Dy5 transformer as a finite element:
i 1 i 2 i 3 i 4 i 5 i 6 i 7 + = G · 2 1 1 a 0 a 0 1 2 1 a a 0 0 1 1 2 0 a a 0 a a 0 a 2 0 0 a 2 0 a a 0 a 2 0 a 2 a 0 a 0 0 a 2 a 2 0 0 0 a 2 a 2 a 2 3 · a 2 · φ 1 φ 2 φ 3 φ 4 φ 5 φ 6 φ 7 + + G · ( 1 ϑ 1 ) · 2 1 1 a 0 a 0 1 2 1 a a 0 0 1 1 2 0 a a 0 a a 0 a 2 0 0 a 2 0 a a 0 a 2 0 a 2 a 0 a 0 0 a 2 a 2 0 0 0 a 2 a 2 a 2 3 · a 2 · φ 1 φ 2 φ 3 φ 4 φ 5 φ 6 φ 7 + H · i 1 i 2 i 3 i 4 i 5 i 6 i 7
Any kind of three-phase transformer group can be modeled by using this procedure.

3.4. Loads Modeled as Finite Elements

Loads in the electric power system are generally defined by the real and reactive power. With the purpose of modeling loads in the electric power system, in this paper, we will derive local systems of equations for a single-phase load and a three-phase star-connected load.
Single-phase load as a finite element is depicted in Figure 8. The single-phase load was modeled as a serial connection of a resistor, capacitor and a coil—in such a way enabling accurate modeling of different types of loads regarding its reactive character.
The starting point for modeling a single-phase load as a finite element are the following voltage differential equations. They are written using the branch current flowing through the single-phase load— i ( t ) (Figure 8):
φ 1 ( t ) φ 2 ( t ) = R · i ( t ) + L · d i ( t ) d t + U c ( t )
i ( t ) = C · d U c ( t ) d t
In the following step, time-integration of these equations is needed, since in that way, instead of ordinary differential equations, we obtain a system of algebraic equations. Again, by performing the numerical time-integration of Equations (57) and (58), from some starting instant t to some ending instant t + , we obtain:
t t + φ 1 ( t ) φ 2 ( t ) R · i ( t ) L · d i ( t ) d t U c ( t ) · d t = 0
t t + i ( t ) C · d U c ( t ) d t · d t = 0
Applying the ϑ -method on Equations (59) and (60), the following is obtained:
φ 1 · ( 1 ϑ ) · Δ t + φ 1 + · ϑ · Δ t φ 2 · ( 1 ϑ ) · Δ t φ 2 + · ϑ · Δ t R · ( 1 ϑ ) · Δ t · i R · ϑ · Δ t · i + L · i + + L · i U c · ( 1 ϑ ) · Δ t U c + · ϑ · Δ t = 0
U c + = U c + ( 1 ϑ ) · Δ t C · i + ϑ · Δ t C · i +
Introducing (62) into (61), we obtain the following after a certain amount of rearrangement:
i + = K · φ 1 + K · φ 2 + + K · ( 1 ϑ 1 ) · φ 1 K · ( 1 ϑ 1 ) · φ 2 K ϑ · U c + N · i
K = ϑ · Δ t L + R · ϑ · Δ t + ϑ 2 · Δ t 2 C
N = L R · ( 1 ϑ ) · Δ t ϑ · ( 1 ϑ ) · Δ t 2 C L + R · ϑ · Δ t + ϑ 2 · Δ t 2 C
Based on the previous equation, we can write the complete local system of equations for the single-phase load using the local nodal currents i 1 ( t ) and i 2 ( t ) :
i 1 i 2 + = K · 1 1 1 1 · φ 1 φ 2 + + K · ( 1 ϑ 1 ) · 1 1 1 1 · φ 1 φ 2 K ϑ 1 1 · U c + N 0 0 N · i 1 i 2
Let us extend the single-phase load to a star-connected three-phase load depicted in Figure 9.
Observing Figure 9 and Equation (63), we can write the following equations for the three-phase load:
i 1 + = K 1 · φ 1 + K 1 · φ 4 + + K 1 · ( 1 ϑ 1 ) · φ 1 K 1 · ( 1 ϑ 1 ) · φ 4 K 1 ϑ · U c 1 + N 1 · i 1
i 2 + = K 2 · φ 2 + K 2 · φ 4 + + K 2 · ( 1 ϑ 1 ) · φ 2 K 2 · ( 1 ϑ 1 ) · φ 4 K 2 ϑ · U c 2 + N 2 · i 2
i 3 + = K 3 · φ 3 + K 3 · φ 4 + + K 3 · ( 1 ϑ 1 ) · φ 3 K 3 · ( 1 ϑ 1 ) · φ 4 K 3 ϑ · U c 3 + N 3 · i 3
i 4 + = ( i 1 + + i 2 + + i 3 + )
In Equations (67)–(69), quantities K i and N i , i = 1 , 2 , 3 are obtained from Equations (64) and (65) by introducing appropriate values of resistance, inductance and capacitance for that particular branch.
Previous system of equations can be written in matrix form which is actually necessary to conform to the finite element method terminology:
i 1 i 2 i 3 i 4 + = K 1 0 0 K 1 0 K 2 0 K 2 0 0 K 3 K 3 K 1 K 2 K 3 K 1 + K 2 + K 3 · φ 1 φ 2 φ 3 φ 4 + + ( 1 ϑ 1 ) · K 1 0 0 K 1 0 K 2 0 K 2 0 0 K 3 K 3 K 1 K 2 K 3 K 1 + K 2 + K 3 · φ 1 φ 2 φ 3 φ 4 1 ϑ · K 1 0 0 0 K 2 0 0 0 K 3 K 1 K 2 K 3 · U c 1 U c 2 U c 3 + N 1 0 0 0 N 2 0 0 0 N 3 N 1 N 2 N 3 · i 1 i 2 i 3
According to Equation (62), expressions for voltages on each of the capacitors of the three-phase star-connected load can now be written as
U c 1 + = U c 1 + ( 1 ϑ ) · Δ t C 1 · i 1 + ϑ · Δ t C 1 · i 1 +
U c 2 + = U c 2 + ( 1 ϑ ) · Δ t C 2 · i 2 + ϑ · Δ t C 2 · i 2 +
U c 3 + = U c 3 + ( 1 ϑ ) · Δ t C 3 · i 3 + ϑ · Δ t C 3 · i 3 +
The system of equations given by (71) represents the complete local system of equations of the finite element of the three-phase star-connected load. Equations (72)–(74) are used for computing the capacitor voltages at the end of the time interval.

4. Test Example

Based on the presented theory, a MATLAB program for power system electromagnetic transient analysis has been developed. To verify the proposed method, electromagnetic transient analysis of the observed power system, shown in Figure 10, has been analyzed. The observed power system, shown in Figure 10, consists of four generators connected via four block-transformers (Dy5) to the power system, four three-phase transmission lines, load and a network equivalent.
  • The generator data are as follows:
  • S n = 120 MVA, P n = 108 MW, Q n = 52.3 MVAr;
  • U n = 14.4 kV, I n = 4811 A;
  • I f 0 = 690 A;
  • I f n = 1192 A.
  • The available parameters of generators are as follows:
  • x d = 1.01 p.u., x q = 0.666 p.u., x q = 0.287 p.u., x d = 0.421 p.u., x d = 0.268 p.u.;
  • x σ = 0.198 p.u., r = 0.00236 p.u., T d = 3.4 s, T d 0 = 7.5 s, T d = 0.0554 s; T q = 0.0876 s, T m G = 8.63 s.
  • Turbine governor data are as follows:
  • T p = 0.02 s, T r = 4.5 s, T w = 2.9 s, T g = 0.5 s;
  • q v _ o p e n = 0.15 , q v _ c l o s e = 0.08 , q m a x = 0.8 , q m i n = 0.1 ;
  • σ = 0.04 , δ 1 = 0.378 , k p = 15 , k i = 0.8 , k d = 1.5 , T d = 5 s;
  • a 11 = 0.5 , a 13 = 1 , a 21 = 1.5 , a 23 = 1 .
  • Excitation system data are as follows:
  • T r = 0.02 s, T a = 0.001 s, T e = 0.1 s, T f = 0.1 s;
  • k r = 1 , k a = 50 , k e = 1 , k f = 0.001 .
  • Step-up transformer data are as follows:
  • S n = 120 MVA, D y 5 ;
  • U n 1 / U n 2 = 10.5 / 220 kV, u k = 12.06 % , P k = 270 kW.
  • Load data are as follows:
  • P n = 120 MW, Q n = 0 , U n = 220 kV.
  • Network equivalent data are as follows:
  • U n = 220 kV, I k 3 = 15.15 kA, I k 1 = 16.46 kA.
Three phase transmission line data for a, b and c phases are as follows: = 41 km
R = 0.18 0.09 0.09 0.09 0.18 0.09 0.09 0.09 0.18 ( Ω / km ) L = 2.2 0.796 0.796 0.796 2.2 0.796 0.796 0.796 2.2 ( mH / km ) C = 7.76 0.76 0.76 0.76 7.76 0.76 0.76 0.76 7.76 ( nF / km ) .
In order to perform the FET procedure, the observed power system (Figure 10) has been divided into 15 finite elements (Figure 11). According to FET terminology, each power system element (generator, transformer, etc.) is represented as a finite element where the local matrices and vectors suitable for assembling the procedure with the other power system elements have been defined by corresponding local systems of equations.
To perform the assembling procedure, based on Figure 11, a connectivity matrix with 27 global nodes can be defined (Table 1).
Based on the connectivity matrix, the global system of equations is obtained by assembling the local system of equations of each finite element according to the standard assembling procedure [10,11]. In the first step, initial conditions are initialized based on operating modes of each generator. At the beginning of the simulation, generator 1 and generator 4 are in the operating mode with P = 108 MW ( 0.9 p . u . ) and Q = 0 MVAr ( 0 p . u . ) , while generator 2 and generator 3 are in operating mode with P = 84 MW ( 0.7 p . u . ) and Q = 52.3 MVAr ( 0.436 p . u . ) . The initialization of active power (P) is accomplished using the governing system where the referent value of the generator mechanical moment ( m m ) is defined. Initialization of reactive power (Q) is accomplished using the excitation system where referent field voltage ( u f ) is defined. At the moment t = 0.8 s, a three-phase fault is set to occur at the middle of transmission line 4 (Figure 10). A three-phase fault has been set by imposing boundary conditions at t = 0.8 s (potentials of global nodes 25, 26 and 27 are set to zero). The duration of the fault is taken to be 100 ms.
In order to demonstrate the accuracy of the proposed method, results obtained by the proposed numerical method have been compared to the solutions obtained by the well-known and established EMTP simulation tool where very good correlation has been demonstrated. Figure 12 and Figure 13 show the currents of phases a, b and c in the transmission line 4 from the direction of busbars A and B in the case of a temporary three-phase short circuit calculated using the FET procedure and compared with the EMTP software package.
As can be seen from Figure 12 and Figure 13, subtransient short-circuit currents of each phase depend on the moment of occurrence of the three-phase short circuit, that is, on the value of the voltage of that phase at the moment of occurrence of the three-phase short circuit. It is a well-known theoretical fact that the largest subtransient current in each phase will appear at the moment when the voltage of that phase is zero. It can be seen that the highest value of the subtransient current occurs in phase b regardless of whether the current has direction from busbar A or B. We can conclude that the voltage of phase b was the lowest at the moment when the three-phase short circuit occurred. In addition, a drop in the direct current (DC) component of the three-phase short-circuit current in all three phases is also clearly visible.
Figure 14 shows the voltage of phase a at the busbar A in the case of a temporary three-phase short circuit calculated using the FET procedure and compared with the EMTP software package.
As can be seen from Figure 14, the voltage of phase a on busbar A dropped to 30 % of its value due to the occurrence of a three-phase short circuit. Note here that voltage drops occurred on busbars B, C, and D; although, they were not shown in this example. The smallest voltage drop is naturally on busbar C since it is connected directly to a relatively strong equivalent network.
Figure 15 shows the field voltage ( u f ) of generator 1 obtained using the excitation system (Appendix B) and compared with the EMTP software package.
Figure 16 shows the rotor angle ( δ ) of all generators obtained using Equation (16) of the generator finite element and compared with the EMTP software package.

5. Discussion

After an exhaustive theoretical derivation of the proposed finite-element-based numerical model for the analysis of electromagnetic transients in the electric power system, a scenario was chosen to demonstrate the calculation possibilities of the proposed model. In the proposed scenario, an electromagnetic transient analysis is performed of the considered power system which consists of four generators connected via four block-transformers (Dy5) to the power system, four three-phase transmission lines, load and network equivalents. At the moment t = 0.8 s , a three-phase fault of duration 100 ms has been simulated in the middle of transmission line 4.
It is a well-known fact that for a relatively short (100 ms) three-phase short circuit, the turbine regulator plays no role in this electromagnetic transient phenomenon due to its extremely large time constants. Therefore, the previously defined mechanical moment ( m m ) was considered constant for all generators during the entire transient event. As can be seen from the attached Figure 12 and Figure 13, the current values of the three-phase short circuit are significantly larger compared to the rated current before the occurrence of the short circuit, regardless of whether the short circuit location is flowing from busbar A or B. As shown in Figure 14, during the short circuit duration, the voltage at busbar A dropped to cca 30 % of its rated value. This significantly affected the rotor angle of generators 1 and 2, as can be seen in Figure 16a,c. Although we did not show them in the Figures, the fact is that there was a voltage drop on busbars B, C and D. Voltage drop on busbar B (see Figure 10) is significantly affected by the rotor angles of generators 3 and 4 as shown in Figure 16b,d. Figure 15 shows the excitation regulator response for generator 1 during the transient event. The fact is that due to its fast time constants, the excitation regulator significantly contributes to maintaining the stability of generator 1, unlike the slow response of the turbine regulator. The same applies to the other three generators. It is also evident that generators 1 and 4 have different load angle responses, given that they have identical parameters as well as initial active and reactive power. The reason is that generator 4 is much closer to the equivalent grid, resulting in less load angle slip due to the smaller voltage drop at busbar B during the first swing. The same applies to generator 2 compared to generator 1. In the case of a longer lasting transient three-phase short circuit, the situation would certainly be more unfavorable, requiring much more time to establish a new steady state. In perspective, we could also consider the action of protection by simulating the disconnection of line 4 from the system. The most unfavorable situation from the perspective of generators 1 and 2 is the three-phase short circuit at busbar A, where they lose interconnection with the rest of the grid. The same applies to generators 3 and 4 in the case of a three-phase short circuit at busbar B.

6. Conclusions

A numerical model for the analysis of electromagnetic transients in a power system, based on the finite element technique, was developed. In this paper, we modeled only the basic elements of the power system. These include numerical models of the synchronous generator (including excitation system and the turbine governing system), transformer, network equivalent, transmission line and load. It is important to emphasize that it is possible to obtain a numerical model of any element of the power system. It is only important to know its mathematical model to define it as a finite element with a certain number of local nodes using the same procedure as outlined in this paper. This new approach enables, if necessary, more detailed modeling of a particular part of the power grid, so that depending on the type of transient phenomena, different mathematical and therefore numerical models of a particular part of the power system can be defined. Based on the presented theory, a MATLAB program for electromagnetic transient analysis of the power system has been developed. The accuracy of the computer program has been validated by comparing it to results obtained by a well-known and established EMTP simulation tool. Excellent agreement of results was found thus confirming the accuracy of the proposed model. EMTP was used as a benchmark tool since it is widely used in scientific and engineering communities and considered to be a most accurate and robust simulation tool.

Author Contributions

Conceptualization, I.J.-G. and D.L.; methodology, I.J.-G.; software, I.J.-G., D.L. and I.K.; validation, I.J.-G. and D.L.; formal analysis, I.J.-G.; investigation, I.J.-G.; writing—original draft preparation, I.J.-G. and D.L.; writing—review and editing, D.L. and I.K.; visualization, I.K. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

Nomenclature

The following quantities and their abbreviations are used in this manuscript:
a , b , c indexes associated with the armature winding on the stator
a 11 , a 13 , a 21 turbine coeficients
a 23 turbine gain
atransformation ratio of single-phase transformer
Csingle-phase load capacitance
D , Q indexes associated with the damping winding on the rotor
findex associated with the excitation winding on the rotor
I f 0 open circuit generator field coil current
I f n nominal generator field coil current
I k 3 three-phase subtransient short circuit current
I k 1 single-phase subtransient short circuit current
i a phase a generator current
i b phase b generator current
i c phase c generator current
i f field coil current
i D damper winding current in axes d
i Q damper winding current in axes q
i p ( t ) current through the primary winding
i s ( t ) current through the secondary winding reduced to the primary side
k p PID controller proportional gain
k i PID controller integral gain
k d PID controller derivative gain
k a , k r , k f , k e regulator amplification factors
three-phase transmission line length
L a ( γ ) phase a stator inductance depending on rotor position
L b ( γ ) phase b stator inductance depending on rotor position
L c ( γ ) phase c stator inductance depending on rotor position
L f ( γ ) field coil inductance depending on angle of rotor position
L D ( γ ) damper winding inductance in axis d depending on angle of rotor position
L Q ( γ ) damper winding inductance in axis q depending on angle of rotor position
L i j ( γ ) mutual inductances between windings i and j depending on rotor position (i, j) = a, b, c, f, D, Q
Lsingle-phase load inductance
L m equivalent network mutual inductance in the natural a, b, c system
L v equivalent network inductance in the natural a, b, c system
L σ 1 leakage inductance of the primary winding
L σ 2 leakage inductance of the secondary winding reduced to the primary side
m m mechanical moment of the generator in p.u.
m e l electrical moment of the generator in p.u.
P k step up transformer short-circuit active power
P n nominal active power
P r e f referent value of mechanical power
P m mechanical power
Q n nominal reactive power
q v _ o p e n gate max opening speed
q v _ c l o s e gate min closing speed
q m a x maximal gate position
q m i n minimal gate position
rgenerator armature (stator) resistance
r a phase a stator resistance
r b phase b stator resistance
r c phase c stator resistance
r f field coil resistance
r D damper winding resistance in axis d
r Q damper winding resistance in axis q
Rsingle-phase load resistance
R 1 resistance of the primary winding
R 2 resistance of the secondary winding reduced to the primary side
sthe synchronous generator rotor slip in p.u.
S k 3 three-phase subtransient short circuit power
S k 1 single-phase subtransient short circuit power
S n nominal apparent power
ttime
T a , T r , T f , T e time constants
T d PID controller derivative time constant
T m the generator mechanical time constant in rad
T p pilot time constant
T r transient droop time constant
T w water inertia time constant
T g gate time constant
T m G the generator mechanical time constant
T d , T d 0 , T d , T q transient and subtransient open- circuit time constants in d and q axes
u a phase a generator voltage
u b phase b generator voltage
u c phase c generator voltage
u f field coil voltage
u f 0 initial field voltage
U n nominal voltage
U c ( t ) single-phase load capacitance voltage
U p ( t ) voltage on the primary winding of the single-phase transformer
U s ( t ) voltage on the secondary winding of the single-phase transformer
U s ( t ) voltage on the secondary winding of the single-phase transformer reduced to the primary side
U n 1 / U n 2 step up transformer nominal voltage ratio
u k step up transformer short-circuit voltage in in percentages
v T , v f , v R , v C regulator signals
v t phase voltage
v R E F referent phase voltage
x d , x q , x d , x d , x q synchronous, transient and subtransient reactances in d and q axes
x σ synchronous generator leakage reactance
X m d direct sequence network reactance
X m i inverse sequence network reactance
X m 0 zero sequence network reactance
X v , X a , X b , X c equivalent network reactance’s in the natural a, b, c system
X m , X a b , X a c , X b c equivalent network mutual reactance’s in the natural a, b, c system
A g l o b global matrix formed by assembling all local matrices using the connectivity matrix
A l o c local matrix of the finite element
A inductance matrix
C capacitance matrix
R resistance matrix
B g l o b global vector formed by assembling all local vectors using the connectivity matrix
B l o c local vector of the finite element
i l o c + vector of local nodal currents at the end of the integration interval
i 1 vector of generator stator winding currents
i 2 vector of generator rotor winding currents
i vector of currents
u 1 vector of generator stator winding voltages
u 2 vector of generator rotor winding voltages
u vector of voltages
γ electrical angle of rotor position
δ 1 temporary droop coefficient
δ the rotor angle in rad
Δ t time increment
ϑ interpolation factor
σ permanent droop coefficient
φ vector of potentials
φ g l o b + vector of global potentials at the end of the integration interval
φ l o c + vector of local nodal potentials at the end of the integration interval
Ψ a phase a stator flux linkage
Ψ b phase b stator flux linkage
Ψ c phase c stator flux linkage
Ψ f field coil flux linkage
Ψ D damper winding flux linkage in axis d
Ψ Q damper winding flux linkage in axis q
Ψ magnetic flux matrix
ω r rotor angular frequency in p.u.

Appendix A. Model of Turbine Governor

The time-domain numerical model of a turbine governor has been developed in [15] and incorporated into the originally developed synchronous generator finite element model. According to [15], the block scheme of the employed turbine control system is given in Figure A1.
Figure A1. Block scheme of a turbine control system.
Figure A1. Block scheme of a turbine control system.
Energies 17 02517 g0a1
The time-domain numerical model of a turbine governor has been derived in detail in [15]. As in this paper, differential equations were reduced to algebraic equations by utilizing the ϑ -method. The resulting equations are given here in a more compact form, and along with the block scheme given in Figure A1, they form a complete system of algebraic equations of the turbine governor used in the finite element model of the synchronous generator. The meaning of used variables are given in the nomenclature section for clarity purposes.
P s c + = T d ( 1 ϑ ) · Δ t · P s c + k d · P s + ( 1 ϑ ) · Δ t · P s T d + ϑ · Δ t
q v + = T g · T p T g · ( 1 ϑ ) · Δ t · q v + ϑ · Δ t · P s 1 + + ( 1 ϑ ) · Δ t · P s 1 T g · T p + T g · ϑ · Δ t
q + = q + ϑ · Δ t · q v + + ( 1 ϑ ) · Δ t · q v
P s b + = P s b + k i · ϑ · Δ t · P s + + k i f · ( 1 ϑ ) · Δ t · P s
P g σ 2 + = T r ( 1 ϑ ) · Δ t · P g σ 2 + δ 1 · T r · q + + δ 1 · T r · q T r + ϑ · Δ t
P m + = a 11 · T w ( 1 ϑ ) · Δ t · P m + a 23 · ϑ · Δ t + a 23 · a 11 · T w a 13 · a 21 · T w · q + ϑ · Δ t + a 11 · T w + a 23 · ( 1 ϑ ) · Δ t + a 13 · a 21 · T w a 23 · a 11 · T w · q ϑ · Δ t + a 11 · T w

Appendix B. Model of the Excitation System

As was the case with the turbine governor, the time-domain numerical model of synchronous generator excitation system has been developed and incorporated into synchronous generator finite element model. According to [24], the block scheme of excitation system model is given in Figure A2.
Figure A2. Block scheme of synchronous generator excitation system.
Figure A2. Block scheme of synchronous generator excitation system.
Energies 17 02517 g0a2
The time-domain numerical model of the excitation system has been derived in detail in [15]. The resulting algebraic equations are given here in a more compact form, and along with the block scheme given in Figure A2, they form a complete system of algebraic equations of the excitation system used in the finite element model of the synchronous generator. The meaning of used variables are again given in the nomenclature section for clarity purposes.
v T + = T r ( 1 ϑ ) · Δ t · v T + k r · ϑ · Δ t · v t + + k r · ( 1 ϑ ) · Δ t · v t T r + ϑ · Δ t
v R + = T a ( 1 ϑ ) · Δ t · v R + k a · ϑ · Δ t · v C + + k a · ( 1 ϑ ) · Δ t · v C T a + ϑ · Δ t
u f + = T e k e · ( 1 ϑ ) · Δ t · u f + ϑ · Δ t · v R + + ( 1 ϑ ) · Δ t · v R T e + k e · ϑ · Δ t
v f + = T f ( 1 ϑ ) · Δ t · v f + k f · u f + k f · u f T f + ϑ · Δ t

References

  1. Anderson, P.M.; Fouad, A.A. Power System Control and Stability, 2nd ed.; IEEE Press: Piscataway, NJ, USA, 2003. [Google Scholar]
  2. Arrillaga, J.; Arnold, C.P.; Harker, B.J. Computer Modelling of Electrical Power Systems; John Wiley and Sons: New York, NY, USA, 1983. [Google Scholar]
  3. Arrillaga, J.; Arnold, C.P. Computer Analysis of Power Systems; John Wiley and Sons: New York, NY, USA, 1990. [Google Scholar]
  4. Dommel, H.W. Digital computer solution of electromagnetic transients in single-end multiphase networks. IEEE Trans. Power Appar. Syst. 1969, 88, 388–389. [Google Scholar] [CrossRef]
  5. Dommel, H.W. Electromagnetic Transients Program Reference Manual (EMTP Theory Book); BPA: Portland, OR, USA, 1986. [Google Scholar]
  6. Woodford, D.A.; Gole, A.M.; Menzies, R.Z. Digital Simulation of DC Links and AC Machines. IEEE Trans. Power Appar. Syst. 1983, 102, 1616–1623. [Google Scholar] [CrossRef]
  7. Nagel, L.W. SPICE2 A Computer Program to Simulate Semiconductor Circuits. Memorandum ERL M520 1975. [Google Scholar]
  8. Mahseredjian, J.; Alvarado, F. Creating an Electromagnetic Transients Program in MATLAB: MatEMTP. IEEE Trans. Power Deliv. 1997, 12, 380–388. [Google Scholar] [CrossRef]
  9. Mahseredjian, J.; Alvarado, F.; Rogers, G.; Long, B. MATLAB’s Power for Power Systems. IEEE J. Comput. Appl. Power 2001, 14, 13–19. [Google Scholar]
  10. Zienkiewicz, O.C.; Taylor, R.L. The Finite Element Method Vol. 1; McGraw-Hill: New York, NY, USA, 1989. [Google Scholar]
  11. Zienkiewicz, O.C.; Taylor, R.L. The Finite Element Method Vol. 2; McGraw-Hill: New York, NY, USA, 1991. [Google Scholar]
  12. Jurić-Grgić, I.; Lucić, R.; Kurtović, M. Power system fault studies using finite element technique. Int. Rev. Model. Simul. 2009, 2, 12–17. [Google Scholar]
  13. Lee, S.Y.; Konrad, A. Lossy Transmission Line Transient Analysis by the Finite Element Method. IEEE Trans. Magn. 1993, 29, 1730–1732. [Google Scholar] [CrossRef]
  14. Lucić, R.; Jurić-Grgić, I.; Kurtović, M. Time domain finite element method analysis of multiconductor transmission lines. Eur. Trans. Electr. Power 2010, 20, 822–832. [Google Scholar] [CrossRef]
  15. Jurić-Grgić, I.; Grulović-Plavljanić, N.; Dabro, M. An Analysis of Power System Transient Stability Using Finite Element Technique. Int. Trans. Electr. Energy Syst. 2019, 29, e2647. [Google Scholar] [CrossRef]
  16. Zienkiewicz, O.C.; Morgan, K. Finite Elements and Approximation; John Wiley & Sons: Hoboken, NJ, USA, 1982. [Google Scholar]
  17. Henschel, S. Analysis of Electromagnetic und Electromechanical Power System Transients with Dynamic Phasors. Ph.D. Dissertation, University of British Colombia, Vancouver, BC, Canada, 1999. [Google Scholar]
  18. Venkatasubramanian, V. Tools for dynamic analysis of the general large power system using time-varying phasors. Int. J. Electr. Power Energy Syst. 1994, 16, 365–376. [Google Scholar] [CrossRef]
  19. Adapa, R.; Reeve, J. A new approach to dynamic analysis of AC networks incorporating detailed modeling of DC systems. II. Application to interaction of DC and weak AC systems. IEEE Trans. Power Deliv. 1988, 3, 2012–2019. [Google Scholar] [CrossRef]
  20. Zambroni de Souza, A.C.; Lopes, B.I.L. Unified computational tool for transient and long-term stability studies. IET Gener. Transm. Distrib. 2009, 3, 173–181. [Google Scholar] [CrossRef]
  21. Zarate, M.R.; Van, C.T.; Federico, M.; Conejo, A.J. Securing transient stability using time-domain simulations within an optimal power flow. IEEE Trans. Power Syst. 2010, 25, 243–253. [Google Scholar] [CrossRef]
  22. Bagde, B.Y.; Umre, B.S.; Dhenuvakonda, K.R. An efficient transient stability-constrained optimal power flow using biogeography-based algorithm. Int. Trans. Electr. Energy Syst. 2018, 28, e2467. [Google Scholar] [CrossRef]
  23. Odun-Ayo, T.; Crow, M.L. An analysis of power system transient stability using stochastic energy functions. Int. Trans. Electr. Energy Syst. 2013, 23, 151–165. [Google Scholar] [CrossRef]
  24. IEEE Standard 421.5-2016; Recommended Practice for Excitation System Models for Power System Stability Studies. IEEE Standard: Piscataway, NJ, USA, 2016.
Figure 1. Electrical schematic of a three-phase synchronous generator with a salient pole rotor.
Figure 1. Electrical schematic of a three-phase synchronous generator with a salient pole rotor.
Energies 17 02517 g001
Figure 2. Synchronous generator as a finite element.
Figure 2. Synchronous generator as a finite element.
Energies 17 02517 g002
Figure 3. Network equivalent as a finite element.
Figure 3. Network equivalent as a finite element.
Energies 17 02517 g003
Figure 4. Single-phase transformer with reduced secondary quantities.
Figure 4. Single-phase transformer with reduced secondary quantities.
Energies 17 02517 g004
Figure 5. Single-phase transformer as one finite element.
Figure 5. Single-phase transformer as one finite element.
Energies 17 02517 g005
Figure 6. Three-phase transformer Dy5 as a finite element.
Figure 6. Three-phase transformer Dy5 as a finite element.
Energies 17 02517 g006
Figure 7. Depiction of a three-phase transformer separated into three pairs of windings.
Figure 7. Depiction of a three-phase transformer separated into three pairs of windings.
Energies 17 02517 g007
Figure 8. Single-phase load as a finite element.
Figure 8. Single-phase load as a finite element.
Energies 17 02517 g008
Figure 9. Three-phase star-connected load.
Figure 9. Three-phase star-connected load.
Energies 17 02517 g009
Figure 10. Single line diagram of an observed electric power system.
Figure 10. Single line diagram of an observed electric power system.
Energies 17 02517 g010
Figure 11. The observed power system divided into finite elements.
Figure 11. The observed power system divided into finite elements.
Energies 17 02517 g011
Figure 12. Currents in transmission line 4 from the direction of busbar A during transient operation for (a) Phase a, (b) Phase b, (c) Phase c.
Figure 12. Currents in transmission line 4 from the direction of busbar A during transient operation for (a) Phase a, (b) Phase b, (c) Phase c.
Energies 17 02517 g012aEnergies 17 02517 g012b
Figure 13. Currents in transmission line 4 from the direction of busbar B during transient operation for (a) Phase a, (b) Phase b, (c) Phase c.
Figure 13. Currents in transmission line 4 from the direction of busbar B during transient operation for (a) Phase a, (b) Phase b, (c) Phase c.
Energies 17 02517 g013
Figure 14. Voltage of phase a at the busbar A during the transient operation.
Figure 14. Voltage of phase a at the busbar A during the transient operation.
Energies 17 02517 g014
Figure 15. Generator 1 field voltage ( u f ) during the transient operation.
Figure 15. Generator 1 field voltage ( u f ) during the transient operation.
Energies 17 02517 g015
Figure 16. Rotor angles ( δ ) during the transient operation for: (a) Generator 1, (b) Generator 4, (c) Generator 2, (d) Generator 3.
Figure 16. Rotor angles ( δ ) during the transient operation for: (a) Generator 1, (b) Generator 4, (c) Generator 2, (d) Generator 3.
Energies 17 02517 g016
Table 1. Connectivity matrix.
Table 1. Connectivity matrix.
Local Nodes
Type Finite Element 1 2 3 4 5 6
Generator 11123
Generator 22456
Generator 33101112
Generator 44131415
Transformer 15123789
Transformer 26456789
Transformer 37101112161718
Transformer 48131415161718
Line 49789252627
Line 410161718252627
Line 211192021222324
Line 312161718192021
Line 113789222324
Equivalent network14192021
Load15222324
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Jurić-Grgić, I.; Lovrić, D.; Krolo, I. Analysis of Power System Electromagnetic Transients Using the Finite Element Technique. Energies 2024, 17, 2517. https://doi.org/10.3390/en17112517

AMA Style

Jurić-Grgić I, Lovrić D, Krolo I. Analysis of Power System Electromagnetic Transients Using the Finite Element Technique. Energies. 2024; 17(11):2517. https://doi.org/10.3390/en17112517

Chicago/Turabian Style

Jurić-Grgić, Ivica, Dino Lovrić, and Ivan Krolo. 2024. "Analysis of Power System Electromagnetic Transients Using the Finite Element Technique" Energies 17, no. 11: 2517. https://doi.org/10.3390/en17112517

APA Style

Jurić-Grgić, I., Lovrić, D., & Krolo, I. (2024). Analysis of Power System Electromagnetic Transients Using the Finite Element Technique. Energies, 17(11), 2517. https://doi.org/10.3390/en17112517

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop