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 phasor-based 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 time-domain simulation has been of interest [1, 2]. Most popular synchronous machine analysis has been established based on a so called two-reaction theory . The two-reaction 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) .
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 sub-windings 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 coil-by-coil basis. The capacitive effect in modeling of the faulted machine is also considered in  and mutual inductance of faulted machine was precisely calculated by use of winding function approach in . The method requires design data such as number of slots and windings structure. The voltage-behind-reactance (VBR) model has recently been proposed for time domain analysis of faulty machine . 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 phasor-based 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 EMTP-type 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 sub-matrix depends on rotor speed that requires re-factorization 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 real-life 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 re-factorization 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 3-phase model of the machine. The network is then a simple Thevenin equivalent. A schematic of the synchronous machine during a single-phase-to-ground fault in phase ‘a’ is shown in Fig. 1.
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.
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 sub-windings, 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, Vabcs, Iabcs 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, Vdqr, Idqr are vectors for rotor fluxes, voltages and currents, respectively. Rs and Rr are diagonal matrices of stator and rotor resistances. The flux linkage relationships and calculation of time-variant inductance matrices, explicitly available in Appendix I of .
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 Leqabc, eabcs are obtained in appendixes II and III of . Eq. (3) is discretized using implicit trapezoidal rule for interfacing the model with the external power system network. Therefore:
Finally, is expressed in terms of the stator currents forming discretized VBR model of faulty machine:
Where Req, Eeq 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.
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 ihis is calculated as ih = GeqEeqhis.
In pure VBR method, since power systems are usually modeled by nodal equations, machine sub-matrix Geq should be inserted into the overall conductance matrix of the system. Moreover, the current history term needs to be injected into the corresponding node . To be noted is that the equivalent conductance matrix, Geq, is rotor-position –dependant resulting in a time-variant machine sub-matrix which has to be re-factorized 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 Diakoptic-Based 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 decomposed-networks are isolated and can be solved independently.
Fig. 4.Network decomposition by Diakoptic technique
Based on Fig. 4, for every ideal CB that is splitting a boundary bus between 2 adjacent sub-networks, following can be written:
where iex is an exchange current between sub-networks. 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 sub-networks according to the Diakoptic concept [19-22]. The tear-off 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 sub-networks. The terminal bus of the faulty machine is torn apart into two buses. The main sub-network is the system circuit without the faulty machine and the other network is the three-phase model of the faulty machine.
Local equations of the sub-networks can be expressed by vector equations:
where, variable yl can be bus voltages or branch currents. These equations are solved through the Newton-Raphson iterative method:
where, Δyk is a vector of incremental changes in voltage or current in kth iteration. The Jacobian matrix of the decomposed network will be written as:
g1 (y1), g2(y2) are the independent sub-network equations. In addition, C1(n1×m) is the incidence matrix that is defined as:
According to LCS, the nominal independent solutions of the sub-networks 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 X0, can be easily obtained as:
By solving (6) for the stator currents of faulty machine, iabc , and overall network equations without incorporating the faulty machine sub-matrix, 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 X0, according to LCS definitions . Then we have:
where X0 is obtained from (16) and Z denotes the current exchange between sub-networks that is calculated from nominal independent solutions. So the sub-vectors of the solution vector X0 may be individually obtained. Note that the sub-vector X1 is the faulty machine variable and subvector X2 is the network variable. X3 represents current increments in supposed ideal CB. At this stage, each state vector can be calculated independently in every time step without re-factorization of overall network matrix. The flowchart of the proposed method is shown in Fig. 5.
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 phasor-based 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 14-Bus 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.
Fig. 6.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 phase-R to 12.5% of phase-T 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.
Fig. 7.Terminal currents of the faulty machine in phase R
Fig. 8.Terminal currents of the faulty machine in phase S
Fig. 9.Terminal currents of the faulty machine in phase T
4.2 Second case study: IEEE-14Bus
In another study case, IEEE-14Bus 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 .
Fig. 10.Case 2, load flow result of IEEE-14Bus 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 Pentium-4 2-GHz processor and 1GB of RAM. Comparisons of the results for stator currents are depicted in Figs. 11-13.
Fig. 11.Current in phase A-terminal side
Fig. 12.Current in phase B
Fig. 13.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 computing-based 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 2-norm formulation employed in  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  to avoid numerical instability due to small number of turns in one of the sub-windings. 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.
Fig. 14.Relative error of IB 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 real-time performance is then easily achievable.
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 three-phase model of a faulty machine. The approach was implemented in stand-alone 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 phasor-based software packages which are using symmetrical component-based transient stability analysis. Moreover, total equations of the network are solved in parallel with the faulty machine equations without re-factorization 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.
The data for lines, transformers and, generators of case 1 are presented in Tables 1, 2 and 3, respectively.
Table 1.Case 1, Cables and lines data
Table 2.Case 1, Transformers data
Table 3.Case 1, Generators data