This paper presents a new approach for simulating the internal faults of synchronous machines using distributed computing and Large Change Sensitivity (LCS) analysis. LCS analysis caters for a parallel solution of 3phase model of a faulted machine within the symmetrical componentbased model of interconnected network. The proposed method considers dynamic behavior of the faulty machine and connected system and tries to accurately solve the synchronous machine’s internal fault conditions in the system. The proposed method is implemented in standalone FORTRANbased phasor software and the results have been compared with available recordings from real networks and precisely simulated faults by use of the ATP/EMTP as a time domain software package. An encouraging correlation between the simulation results using proposed method, ATP simulation and measurements was observed and reported. The simplified approach also enables engineers to quickly investigate their particular cases with a reasonable precision.
1. Introduction
Stator windings in synchronous machines are insulated from the stator core by winding insulation. The material used in stator winding insulation is made up of mainly organic composite materials and thus have a poor thermal and mechanical property in comparison with copper and iron. During years, aging and stresses on the insulation materials weakens the insulations and, depending on the effectiveness of operation and maintenance, insulation failure might happen resulting in stator earth faults. Whereas detailed simulation of such earth faults is a necessity, there has not been a global effort to facilitate stator earth fault simulations in power system simulation toolkits. Moreover, as it is customary for protection engineers to use simplified analytical methods on their application, there has been a lack of simplified engineering approaches to bridge the complex simulation providing reasonably accurate results. Bearing this in mind, in this paper, a methodology has been proposed to incorporate stator faults in phasorbased dynamic analysis algorithms and a simplified method is developed to incorporate stator earth fault calculation to be used in engineering applications.
From the works has been done so far on the subject, experimental tests and extensive timedomain simulation has been of interest
[1
,
2]
. Most popular synchronous machine analysis has been established based on a so called tworeaction theory
[3]
. The tworeaction theory, while absolutely effective in fast and accurate analysis, is not suitable for stator winding short circuit simulation because the conventional dq0 model is no longer applicable. To overcome this shortcoming, the synchronous machine fault representation is carried out in the Phase Domain (PD)
[4]
.
For instance, in
[5

8]
, partitioning of machine windings is used to calculate the faulty machine inductance matrix in phase domain where a faulted winding is considered to be made up of two subwindings with the effective displacement of the magnetic axis. Alternatively in
[9

11]
, the synchronous machine is assumed to be formed from several electric loops and the inductances are calculated on a coilbycoil basis. The capacitive effect in modeling of the faulted machine is also considered in
[12]
and mutual inductance of faulted machine was precisely calculated by use of winding function approach in
[13]
. The method requires design data such as number of slots and windings structure. The voltagebehindreactance (VBR) model has recently been proposed for time domain analysis of faulty machine
[14]
. The model is developed using a PD representation of partitioned winding.
In brief, for the sake of internal fault modeling of synchronous machines in conjunction with external power system network, two methodologies exist. The first approach deploys static Thevenin equivalent circuit of the power system. This method exhibits limited success for studies which require high precisions such as relay coordination. The second is the detailed
abc
phase coordinate model of power system. In latter approach, storage and CPU simulation time grow rapidly with the size of the power system. Since power system networks are often modeled in
abc
phase in time domain software such as ATP and in 012 sequence coordinates in phasorbased software, the interface between faulted machine model and the external power system network develops further challenges. To achieve an interface between the machine model and external network, several models have been recently proposed
[15
,
16]
. In EMTPtype solution
[17
,
18]
, the VBR model of induction machine formulation represents the stator circuit in
abc
phase coordinates and the rotor circuit in
dq
reference frame. The representation employs the benefits of the PD model for direct interface with the
abc
phase coordinate model of external network. However, the machine conductance submatrix depends on rotor speed that requires refactorization of the entire network conductance matrix at every time step as the rotor speed changes. It increases the computational cost especially in large interconnected networks and thus limits the application in reallife engineering simulations.
In this paper an efficient approach is developed to deploy the VBR model of faulted machine to simulate the internal fault condition without refactorization of the network conductance matrix at every time step. To achieve this goal the paper combines the discretized VBR model and distributed simulation concepts (Diakoptic
[19

22]
) and LCS
[23

26]
. In addition to an encouraging match between the ATP simulation, proposed approach and measured data, the simulation speed is considerably higher. This approach is primarily of use for engineers to simulate their cases with accurate and fast simulation using existing data without the need for time consuming modeling and simulations.
2. Three Phase Model of Faulty Synchronous Machine
A widely used method to model internal fault condition of synchronous machines utilizes the 3phase model of the machine. The network is then a simple Thevenin equivalent. A schematic of the synchronous machine during a singlephasetoground fault in phase ‘a’ is shown in
Fig. 1
.
Schematic of the faulty synchronous machine
In order to simulate an internal fault, it is necessary to divide the faulty windings. The details of this partitioning can be found in
[6

8
]. An example of partitioning is shown in
Fig. 2
.
Partitioning of faulty winding
Angle α defines the angular location of the internal node f, which divides the winding into 2 sections. θ, γ
_{1}
, γ
_{2}
are the angular position of rotor at any instant and displacement of magnetic axis of subwindings, respectively.
The voltage equations in the phase reference frame of synchronous machine may be written as:
where

Ψabcs=[ΨatΨbtΨctΨanΨbnΨcn]T

Vabcs=[VatVbtVctVanVbnVcn]T

Iabcs=[IatIbtIctIanIbnIcn]T
Ψ_{abcs}
,
V_{abcs}
,
I_{abcs}
are vectors corresponding to the stator windings fluxes, voltages and currents; the subscripts t and n refer to terminal and neutral sides of the partitioned windings, respectively. P is the d/dt operator;
Ψ_{dqr}
,
V_{dqr}
,
I_{dqr}
are vectors for rotor fluxes, voltages and currents, respectively. R
_{s}
and R
_{r}
are diagonal matrices of stator and rotor resistances. The flux linkage relationships and calculation of timevariant inductance matrices, explicitly available in Appendix I of
[14]
.
The VBR model develops the PD model which expresses all the rotor and stator circuits using abc phase quantities. The model is established upon using the transformation of phase variable to the rotor reference frame and Park’s transformation. After performing algebraic simplifications, the final equations for stator voltages can be expressed as:
where L
_{eqabc}
, e
_{abcs}
are obtained in appendixes II and III of
[14]
. Eq. (3) is discretized using implicit trapezoidal rule for interfacing the model with the external power system network. Therefore:
where:
Finally,
is expressed in terms of the stator currents forming discretized VBR model of faulty machine:
Where
R_{eq}
,
E_{eq his}
are the equivalent resistance matrix and equivalent voltage history term, respectively. The final equivalent network as a single machine that is connected to the power system can be similar to what is presented in
Fig. 3
.
Interface of machine model with power system
It is important now to obtain the numerical solution for the faulty machine model with the external network. Moving forward, the machine branch Eqs. (6) are replaced by the machine nodal expression as:
where,
is the equivalent conductance matrix of faulty machine and the current history term
i_{his}
is calculated as
i_{h}
=
G_{eq}E_{eqhis}
.
In pure VBR method, since power systems are usually modeled by nodal equations, machine submatrix
G_{eq}
should be inserted into the overall conductance matrix of the system. Moreover, the current history term needs to be injected into the corresponding node
[12]
. To be noted is that the equivalent conductance matrix,
G_{eq}
, is rotorposition –dependant resulting in a timevariant machine submatrix which has to be refactorized at every time step of simulation. This requirement significantly increases the overall computational burden of simulations. To overcome this problem, one idea is to solve the internal fault of synchronous machine together with network and generators dynamics utilizing the LCS and network decomposition. We will discuss this approach in following sections.
3. LCS and DiakopticBased Solution for Interfacing Faulted Machine to Network
Assume an electrical network with n buses. Based on the Diakoptic technique, this network can be split into smaller networks which are connected to each other through supposed ideal circuit breakers (CBs), as shown in
Fig. 4
. If the CBs are open the decomposednetworks are isolated and can be solved independently.
Network decomposition by Diakoptic technique
Based on
Fig. 4
, for every ideal CB that is splitting a boundary bus between 2 adjacent subnetworks, following can be written:
where
i_{ex}
is an exchange current between subnetworks. If F is 0 then the CB is open and buses at either side of the CB have different voltages. But, if F is 1, the CB is closed and either side of the boundary buses have same voltage. This technique is usually implemented for parallel processing of large networks to decrease the computation time. This technique can be used for internal fault analysis of the machine network interfacing challenges.
To do so we tear apart the interconnected network into two subnetworks according to the Diakoptic concept
[19

22]
. The tearoff point is the terminal of the faulty machine. Now we can apply LCS and, based on the LCS definitions
[24

26]
, we need to obtain a nominal independent solution of the subnetworks. The terminal bus of the faulty machine is torn apart into two buses. The main subnetwork is the system circuit without the faulty machine and the other network is the threephase model of the faulty machine.
Local equations of the subnetworks can be expressed by vector equations:
where, variable
y_{l}
can be bus voltages or branch currents. These equations are solved through the NewtonRaphson iterative method:
where, Δ
y^{k}
is a vector of incremental changes in voltage or current in k
^{th}
iteration. The Jacobian matrix of the decomposed network will be written as:
where
g
_{1}
(
y
_{1}
),
g
_{2}
(
y
_{2}
) are the independent subnetwork equations. In addition,
C
_{1}
(
n
_{1}
×
m
) is the incidence matrix that is defined as:
According to LCS, the nominal independent solutions of the subnetworks are obtained when the value of F in (11)(12) is set to zero. Because of the Blocked Bordered Diagonal Form (BBDF) of the Jacobian matrix, the state vector
X
_{0}
, can be easily obtained as:
By solving (6) for the stator currents of faulty machine,
i_{abc}
, and overall network equations without incorporating the faulty machine submatrix, the vector of bus voltages can be obtained as:
Now, when the value of F in (11)(12) is set to unity, the interconnected network with faulty machine is obtained and (10) can be expressed by:
It can be rewritten as:
X
can be calculated by updating the nominal independent solution of X
_{0}
, according to LCS definitions
[27]
. Then we have:
where
X_{0}
is obtained from (16) and
Z
denotes the current exchange between subnetworks that is calculated from nominal independent solutions. So the subvectors of the solution vector
X_{0}
may be individually obtained. Note that the subvector
X_{1}
is the faulty machine variable and subvector X
_{2}
is the network variable.
X_{3}
represents current increments in supposed ideal CB. At this stage, each state vector can be calculated independently in every time step without refactorization of overall network matrix. The flowchart of the proposed method is shown in
Fig. 5
.
Flowchart of the proposed method for simulating the synchronous machine internal faults
4. Simulation Results and Discussions
The proposed procedure has been implemented in phasorbased software in FORTRAN and has been applied to two separate cases. The first case is the Karun Network, south west network of Iran (KNIRI). The second is an IEEE 14Bus test system.
 4.1 First case study: KNIRI network
In this section, the proposed approach has been applied to a real power system.
Fig. 6
shows the single line diagram of the network. The data for lines, transformers and generators are presented in Appendix.
Case 1, structure of KNIRI network
An internal fault has been occurred in stator of one of the generating unit in the network (Gen.7) causing forced outage of the unit. The fault has later being investigated during maintenance procedure to be a two phase to ground fault in 25% of phaseR to 12.5% of phaseT in stator windings. The percentage of turns is measured from the phase terminal to the neutral.
The aforementioned method is deployed in symmetrical component based software to simulate this fault. The recorded values of the installed relays are compared with simulation results in
Figs. 7

9
for each phase. The per unit (pu) base of stator p hase currents is in a common 100 MVA base. Apparently, the proposed method results matches quite well with the recorded values.
Terminal currents of the faulty machine in phase R
Terminal currents of the faulty machine in phase S
Terminal currents of the faulty machine in phase T
 4.2 Second case study: IEEE14Bus
In another study case, IEEE14Bus test system is studied in by use of the proposed algorithm. The system has 14 buses, 5 generator units and 20 branches. Single line diagram and the initial conditions for this standard system are shown in
Fig. 10
. The dynamic data of this network is available in
[28]
.
Case 2, load flow result of IEEE14Bus test system
In this study, the machine at Bus 8 is subjected to an internal single phase to ground short circuit in phase A at 10% of the total number of turns at a time 0.1 s, and then removed at 0.3 s, due to protection relay operation. The percentage of turns is measured from the phase terminal to the neutral.
To compare the numerical accuracy of the proposed method with time domain simulation software, the test system was modeled in Alternate Transients Program (ATP). Simulations were run with small time steps, Δ
t
= 1
μ
s using a personal computer (pc) with a Pentium4 2GHz processor and 1GB of RAM. Comparisons of the results for stator currents are depicted in
Figs. 11

13
.
Current in phase Aterminal side
Current in phase B
Current in phase C
As it can be seen from these figures, there is a good agreement between the results from ATP simulation and proposed algorithm. In addition, the measured CPU time for the case studies shows that the proposed method is about 12 times faster than the ATP solution. Proposed distributed computingbased method is the suitable for studies that require quick engineering simulations such as relay coordination.
 4.3 Numerical accuracy
To further compare the numerical accuracy and robustness of the proposed method, the fault condition is simulated for time steps ranging from 1 to 1000 µs. The relative error of the proposed method, calculated using the 2norm formulation employed in
[29]
for varying time steps, is compared to ATP simulations. The 0.1 µs time step of ATP simulation is selected as an exact reference solution. The highest relative error observed in these figures is 3.5% which is considered to be small enough not to affect the accuracy of the method.
 4.4 Effect of fault location
Similar cases were investigated considering different time steps and different internal fault locations to compare the proposed method with ATP solution. The fault location has been adjusted to be within the recommended limits in
[2]
to avoid numerical instability due to small number of turns in one of the subwindings.
Fig. 14
depicts the relative error for varying internal fault locations between ATP and proposed method. It’s worth mentioning that the proposed method has an adequate stable response for varying fault locations.
Relative error of I_{B} for varying internal fault location in phase A
To be noted is that for a goal of achieving 1% relative error (i.e. acceptable accuracy), the time step must be less than 400µs. If this time step (400µs) is selected for simulation runs, the execution time for proposed method will be approximately 20 times faster than 0.1 µs time step. The realtime performance is then easily achievable.
5. Conclusion
In this paper, a novel approach for transient analysis of power systems under internal stator fault conditions of synchronous machines was presented. The method is based on a distributed analysis using Diakoptic, large change sensitivity and a threephase model of a faulty machine. The approach was implemented in standalone software using FORTRAN. The simulation results were compared with the recorded values in a real machine fault. Moreover, the proposed method was compared against ATP simulations. Both comparisons implied that the proposed approach is producing accurate results for internal fault calculation in machines connected to the network. Major advantages of the proposed method are the simplicity in implementation, full compatibility with phasorbased software packages which are using symmetrical componentbased transient stability analysis. Moreover, total equations of the network are solved in parallel with the faulty machine equations without refactorization of entire network matrix. It can simulate fault currents in any branch of the network which is a useful feature for relays coordination and similar engineering applications.
Appendix
The data for lines, transformers and, generators of case 1 are presented in
Tables 1
,
2
and
3
, respectively.
Case 1, Cables and lines data
Case 1, Cables and lines data
Case 1, Transformers data
Case 1, Transformers data
Case 1, Generators data
BIO
Amangaldi Koochaki was born in 1981; received the B.Sc in electrical engineering from the Faculty of Electrical Engineering, University of Tehran, Iran, in 2003; he received M.Sc and Ph.D from Amirkabir University of Technology, Iran, in 2003 and 2010 respectively. Now he is assistant professor in Aliabad Katoul branch of Islamic Azad University, Iran. His research interests are power system analysis, relay coordination and renewable energies.
Reichmeider P. P.
,
Gross C. A.
,
Querry D.
,
Novosel D.
,
Salon S.
2000
“Internal fault in synchronous machines part I: The machine model
“IEEE Trans. Energy Convers.
15
(4)
376 
379
DOI : 10.1109/60.900496
Reichmeider P. P.
,
Gross C. A.
,
Querry D.
,
Novosel D.
,
Salon S.
2000
“Internal fault in synchronous machines part II: Model performance
IEEE Trans. Energy Convers.
15
(4)
380 
383
DOI : 10.1109/60.900497
Park R. H.
1929
“Tworeaction theory of synchronous machines,”
AIEE Trans.
48
(2)
716 
727
Hemmati S.
,
Shokri kojoori S.
,
Saied S.
,
Lipo T.A.
2013
“Modeling and experimental validation of internal shortcircuit fault in salientpole synchronous machines using numerical gap function including stator and rotor core saturation,”
IET Electr. Power Appl.
7
(5)
391 
399
DOI : 10.1049/ietepa.2012.0342
Muthumuni D.
,
McLaren P. G.
,
Dirks E.
,
Pathirana V.
“A synchronous machine model to analyze internal faults,”
Proc. 2001 IEEE Ind. Appl. Conf.
Chicago, IL
Sep. 30Oct. 4, 2001
1595 
1600
Megahed A. I.
,
Malik O. P.
1999
“Simulation of internal faults in synchronous generators,”
IEEE Trans. Energy Convers.
14
(4)
1306 
1311
DOI : 10.1109/60.815064
Megahed A. I.
,
Malik O.P.
1998
“Synchronous generator internal fault computation and experimental verification,”
IET Gener. Transm. Dist.
145
(5)
604 
610
DOI : 10.1049/ipgtd:19982190
Reichmeider P. P.
,
Querry D.
,
Gross A.
,
Novosel D.
,
Salon S.
2000
“Partitioning of synchronous machine winding for internal fault analysis,”
IEEE Trans. Energy Convers.
15
(4)
372 
375
DOI : 10.1109/60.900495
Wang X.
,
Chen S.
,
Wang W.
,
Sun Y.
,
Xu L.
2002
“A study of armature winding internal faults for turbogenerators,”
IEEE Trans. Ind. Appl.
38
(3)
625 
631
DOI : 10.1109/TIA.2002.1003410
Bi D.
,
Wang X.
,
Wang W.
,
Zhu Z. Q.
,
Howe D.
2005
“Improved transient simulation of salientpole synchronous generators with internal and ground faults in the stator winding,”
IEEE Trans. Energy Convers.
20
(1)
128 
134
DOI : 10.1109/TEC.2004.841509
Wang X. H.
,
Sun Y. G.
,
Ouyang B.
,
Wang W. J.
,
Zhu Z. Q.
,
Howe D.
2002
“Transient behavior of salientpole synchronous machines with internal stator winding faults,”
IET Electr. Power Appl.
149
(2)
143 
151
DOI : 10.1049/ipepa:20020068
Lin X.
,
Tian Q.
,
Gao Y.
,
Liu P.
2007
“Studies on the internal fault simulation of a high voltage cable wound generator,”
IEEE Trans. Energy Convers.
22
(2)
240 
249
DOI : 10.1109/TEC.2006.874206
Tu X.
,
Dessaint L.
,
Kahel M. E.
,
Barry A. O.
2006
“A new model of synchronous machine internal faults based on winding distribution,”
IEEE Trans. Ind. Electron.
53
(6)
1818 
1828
DOI : 10.1109/TIE.2006.885125
Rodriguez D. S.
,
Acha E.
2009
“A synchronous generator internal fault model based on the voltagebehind reactance representation,”
IEEE Trans. Energy Convers.
24
(1)
184 
194
DOI : 10.1109/TEC.2008.2005281
Gao W.
,
Solodovnik E. V.
,
Dougal R. A.
2004
“Symbolically aided model development for an induction machine in virtual test bed,”
IEEE Trans. Energy Convers.
19
(1)
125 
135
DOI : 10.1109/TEC.2003.822316
Wang L.
,
Jatskevich J.
2010
“Approximate voltagebehindreactance induction machine model for efficient interface with EMTP network solution,”
IEEE Trans. Power Syst.
25
(2)
1016 
1031
DOI : 10.1109/TPWRS.2009.2034526
Wang L.
,
Jatskevich J.
2006
“A voltagebehindreactance synchronous machine model for the EMTPtype solution,”
IEEE Trans. Power Syst.
21
(4)
1539 
1549
DOI : 10.1109/TPWRS.2006.883670
Wang L.
,
Jatskevich J.
,
Wang C.
,
Li P.
2008
“A voltagebehindreactance induction machine model for the EMTPtype solution,”
IEEE Trans. Power Syst.
23
(3)
1226 
1238
DOI : 10.1109/TPWRS.2008.926423
Kron G.
1952
“Tensorial analysis of integrated transmission systems, part III – the primitive division,”
AIEE Trans.
71
(1)
814 
821
Happ H.H.
1980
“Piecewise Methods and Applications to Power Systems”
John Wiley & Sons
New York
Rohrer R. A.
1988
“Circuit partitioning simplified,”
IEEE Trans. Circ Syst.
35
(1)
2 
5
DOI : 10.1109/31.1694
Kaypmaz A.
,
Tokad Y.
,
Usta U.
“Diakoptics – a parallel processing approach for the analysis of largescale power systems,”
Proc. IEEE 7th Mediterranean Electro technical Conf.
Antalya, Turkey
April 1214, 1994
1069 
1072
Iordache M.
,
Iordache B.
“Generalized diakoptic analysis for large scale piecewiselinear nonlinear electrical circuits,”
IEEE Int. Symp. on Circuits and Systems
London, England
30 May2 June 1994
41 
44
Stenbakken G. N.
,
Starzyk J. A.
1992
“Diakoptic and large change sensitivity analysis,”
IET Circ. Devices Syst.
139
(1)
114 
118
Kalantari A.
,
Kouhsari S.M.
2005
“A distributed computing approach for power system analysis,”
Iran J Sci Technol B
29
(4)
399 
413
Esmaeili S.
,
Kouhsari S. M.
2007
“A distributed simulation based approach for detailed and decentralized power system transient stability,”
Electr. Power Syst. Res.
77
673 
684
DOI : 10.1016/j.epsr.2006.06.008
Koochaki A.
,
Kouhsari S. M.
2010
“Detailed simulation of transformer internal fault in power system by diacoptical concept,”
Advances in electrical and computer engineering
10
(3)
48 
54
Anderson PM.
,
Fouad AA
1977
Power system control and stability
Iowa State University
Ames, IA
Zhu W.
,
Pekarek S.
,
Jatskevich J.
,
Wasynczuk O.
,
Delisle D.
2005
“A modelintheloop interface to represent source dynamics in a hardware based dc distribution system,”
IEEE Trans. Power Electron.
20
(2)
438 
445
DOI : 10.1109/TPEL.2004.842973