1 Transient Stability Problem

The six steps for computing transient stability problem

1. Obtain the generator data for the generator model.

2. Obtain the pre-fault network state information, e.g. network impedances. (The overdetermined network will be reduced before forming the admittance matrix.)

3. Use loadflow analysis to compute initial data of the generators’ terminal voltages and the mechanical power.

4. Obtain the post-fault network state information, e.g. fault location, fault clearing time and breaker reclosings time.

5. Obtain the timing of the simulated system, e.g. time of switching and the length of simulation time.

6. Compute the rotor angle using DAE solver.

The last three steps are massively parallel for screening fault scenarios. All the data should be normalize to a common base unit. A system base of 100 MVA is generally used. Simulation is started before switching occurs in the network to safeguard the simulation from inconsistent data. In addition to the generator equation, we needed the followings to simulate a complete system

• load or network constraints

• the excitation system

• mechanical torque

Detail model of these variables are nonlinear. However, for short transient stability simulation, a linearized model can sometime be used. A linearize model can be used (with care) to simulate the network’s oscillating behavior. The discontinuity across the network is caused by switching and fault. Limits of controllers contributes to the local discontinuity for a particular device. Reduced network model is used to study the effect of a fault in the local area. Remote network is transformed to equivalent circuit and classical generator model is used for remote generator.

2 Classical Generator Model

Assumptions used for classical model

• The mechanical power input remains constant during the period of the transient. (P_{m} = Constant)

• Damping or asynchronous power is negligible (This can be simulated by adding a damping term).

• The synchronous machine can be represented (electrically) by a constant voltage source behind a transient reactance. (This implies the magnitude of the internal machine voltages remain constant for the simulation time.)

• The mechanical angle of the synchronous machine rotor coincides with the electrical phase angle of the voltage behind transient reactance. (The delta angle of the internal machine voltages induced by the transient reactance is the delta rotor angle.)

• If a local load is fed at the terminal voltage of the machine, it can be represented by a constant impedance (or admittance) to neutral. (No load change except the change cause by the faulty component.)

The change in delivered power is reflected by the change in current flow. This model is designed only for the first swing (ie. accurate up to 1 to 2 sec simulation time).

3 Classical Generator Model for Multimachine System

Assumptions applying classical model to multimachine system

• Mechanical power input is constant.

• Damping or asynchronous power is negligible (This can be simulated by adding a damping term to each simulated generator).

• Constant-voltage-behind-transient-reactance model for the synchronous machines is valid.

• The mechanical rotor angle of a machine coincides with the angle of the voltage behind the transient reactance.

• Loads are represented by passive impedances. (No non-linear load.)

These assumptions only true for studying first swing period.

4 Classical Generator equations for Multimachine System

For i = 1, 2, ⋯, n

\begin{aligned}\frac{2H_{i}}{ω_{R}} \frac{dω_{i}}{dt} + D_{i}ω_{i} & = & P_{mi} − [E_{i}^{2}G_{ii} + ∑_{j=1,j≠i}^{n} E_{i}E_{j} cos(θ_{ij} − δ_{i} + δ_{j})] \\ \frac{dδ_{i}}{dt} & = & ω_{i} − ω_{R} \end{aligned}

At (t = 0^{−})P_{mi0} = P_{ei0}

P_{mi0} = E_{i}^{2}G_{ii0} + ∑_{j=1,j≠i}^{n} E_{i}E_{j}Y_{ij0} cos(θ_{ij0} − δ_{i0} + δ_{j0})

G_{ij} is the real part of the i th, j th entry of the Y matrix. E_{i} is the voltage of the node i. In the simulation, this value will obtain by a load flow calculation. D is the damping factor.

5 2-axis Generator Model

Assumptions used for 2-axis generator model

• all faults are three phase fault.

• neglected any subtransient effects

6 The equation for 2-axis Generator Model in state-space from

[ i̇_{d} ; i̇_{F} ; i̇_{D} ; i̇_{q} ; i̇_{Q} ; ω̇ ; δ̇ ] = [ \begin{matrix}−L^{−1}(R + ωN) & 0 \\ \begin{matrix}−\frac{L_{d}i_{q}}{3τ_{j}} & −\frac{kM_{F} i_{q}}{3τj} & −\frac{kM_{D}i_{q}}{3τ_{j}} & \frac{L_{q}i_{d}}{3τ_{j}} & \frac{kM_{Q}i_{d}}{3τ_{j}} \\ 0 & 0 & 0 & 0 & 0 \end{matrix} & \begin{matrix}−\frac{D}{3τ_{j}} & 0 \\ 1 & 0 \end{matrix} \end{matrix} ] [ \begin{matrix}i_{d} \\ i_{F} \\ i_{D} \\ i_{q} \\ i_{Q} \\ ω \\ δ \end{matrix} ] + [ −L^{−1}v ; \frac{T_{m}}{3τ_{j}} ; −1 ]

L = [ \begin{matrix}L_{d} & kM_{F} & kM_{D} & 0 & 0 \\ kM_{F} & L_{F} & M_{R} & 0 & 0 \\ kM_{D} & M_{R} & L_{D} & 0 & 0 \\ 0 & 0 & 0 & L_{q} & kM_{Q} \\ 0 & 0 & 0 & kM_{Q} & L_{Q} \end{matrix} ]

R + ωN = [ \begin{matrix}r & 0 & 0 & ωL_{q} & ωkM_{q} \\ 0 & r_{F} & 0 & 0 & 0 \\ 0 & 0 & r_{D} & 0 & 0 \\ −ωL_{d} & −ωkM_{F} & −ωkM_{D} & r & 0 \\ 0 & 0 & 0 & 0 & r_{Q} \end{matrix} ]

This set of equation is the most general form of the two-axis model. Simplified version of the above model is often used, when the model is used to simulate a specific problem. In addition to the above system of equation, we needed the followings to simulate a complete system

• load or network constraints

• the excitation system

• mechanical torque

An alternate state-space model, which is called flux linkages model, with the state variables, λ_{d}, λ_{F}, λ_{D}, λ_{q} and λ_{Q}. The relationship between the current and flux linkage is give by

[ \begin{matrix}λ_{d} \\ λ_{F} \\ λ_{D} \\ λ_{q} \\ λ_{Q} \end{matrix} ] = [ \begin{matrix}L_{d} & 0 & kM_{F} & kM_{D} & 0 \\ 0 & L_{q} & 0 & 0 & kM_{Q} \\ kM_{F} & 0 & L_{F} & M_{R} & 0 \\ kM_{D} & 0 & M_{R} & L_{D} & 0 \\ 0 & kM_{Q} & 0 & 0 & L_{Q} \end{matrix} ] [ \begin{matrix}i_{d} \\ i_{F} \\ i_{D} \\ i_{q} \\ i_{Q} \end{matrix} ]

The system of equations is for generic generators. Consequently, some variables may not appropriately used to describe a specific type of generators. Those is generally set to zero.

Figure 1: 2-axis machine model 4

7 EPRI’s ETMSP

ETMSP uses explicit RK4 at time of switching with step size zero to handle the discontinuity. Then ETMSP can switch over to implicit trapezoidal rule. ETMSP ignores the local discontinuity (This reduce the accuracy of the simulation). The assumptions of generator model used in ETMSP

• Up to three rotor windings at each axis can be represented

• d and q axes saturation can be separately represented

• Subtransient saliency is fully represented

• An interface is provided for excitation system and speed governing system models

• Rate of change of stator flux linkages and the effect of speed variations are neglected in the stator voltage equations (Complexity of the model is reduced).

Schematic of the 2 axis generator model used by ETMSP is shown in figure 2.

7.1 State Equation

\frac{dΨ_{f}}{dt} = −ω_{0}r_{f}i_{f} + ω_{0}V_{f}
\frac{dΨ_{kd1}}{dt} = −ω_{0}r_{kd1}i_{kd1}
\frac{dΨ_{kd2}}{dt} = −ω_{0}r_{kd2}i_{kd2}
\frac{dΨ_{kq1}}{dt} = −ω_{0}r_{kq1}i_{kq1}
\frac{dΨ_{kq2}}{dt} = −ω_{0}r_{kq2}i_{kq2}
\frac{dΨ_{kq3}}{dt} = −ω_{0}r_{kq3}i_{kq3}
\frac{dω_{m}}{dt} = \frac{ω_{0}}{2H} [T_{m} − (Ψ_{d}i_{q} − Ψ_{q}i_{d})]
\frac{dδ}{dt}  = ω_{m}

t - time (sec) δ - rotor angle (rad) ω_{m} - rotor speed (rad/sec) V_{f} - the field voltage determined by the excitation control system (p.u.) V_{f} is constant, if no excitation control T_{m} - the mechanical torque determined by the speed governing system. (p.u.) T_{m}is constant, if no speed governor

ω_{m}, δ, Ψ_{f}, Ψ_{kd1}, Ψ_{kd2}, Ψ_{kq1}, Ψ_{kq2}, Ψ_{kq3} are the state variable of the generator model. The damper winding may not present.

7.2 Stator Equations

\begin{aligned}V_{d} & = & −r_{a}i_{d} − Ψ_{q} \\ V_{q} & = & −r_{a}i_{q} − Ψ_{d} \end{aligned}

V_{d} - machine terminal voltage in machine d reference frame. V_{q} - machine terminal voltage in machine q reference frame. i_{d} - machine terminal current in machine d reference frame. i_{q} - machine terminal current in machine q reference frame.

Figure 2: EPRI detail synchronous machine model 6

7.3 Flux Linkage Equation

[ \begin{matrix}Ψ_{d} \\ Ψ_{f} \\ Ψ_{kd1} \\ Ψ_{kd2} \end{matrix} ] = [ \begin{matrix}X_{d} & X_{ad} & X_{ad} & X_{ad} \\ X_{ad} & X_{f} & X_{fad1} & X_{fad2} \\ X_{ad} & X_{fad1} & X_{ad1} & X_{ad12} \\ X_{ad} & X_{fad2} & X_{ad12} & X_{ad2} \end{matrix} ] [ \begin{matrix}−i_{d} \\ i_{f} \\ i_{kd1} \\ i_{kd2} \end{matrix} ]

[ \begin{matrix}Ψ_{q} \\ Ψ_{kq1} \\ Ψ_{kq2} \\ Ψ_{kq2} \end{matrix} ] = [ \begin{matrix}X_{q} & X_{aq} & X_{aq} & X_{aq} \\ X_{aq} & X_{aq1} & X_{aq} & X_{aq} \\ X_{aq} & X_{aq} & X_{aq2} & X_{aq} \\ X_{aq} & X_{aq} & X_{aq} & X_{aq3} \end{matrix} ] [ \begin{matrix}−i_{q} \\ i_{kq1} \\ i_{kq2} \\ i_{kq3} \end{matrix} ]

\begin{aligned}X_{d} & = & X_{ad} + X_{l} \\ X_{f} & = & X_{ad} + X_{fd} + X_{c1} + X_{c2} \\ X_{fad1} & = & X_{ad} + X_{c1} \\ X_{fad2} & = & X_{ad} + X_{c1} + X_{c2} \\ X_{ad1} & = & X_{ad} + X_{kdl} + X_{c1} \\ X_{ad12} & = & X_{ad} + X_{c1} \\ X_{ad2} & = & X_{ad} + Z_{kd2} + X_{c1} + X_{c2} \\ X_{q} & = & X_{aq} + X_{l} \\ X_{aq1} & = & X_{aq} + X_{kq1} \\ X_{aq2} & = & X_{aq} + X_{kq2} \\ X_{aq3} & = & X_{aq} + X_{kq3} \end{aligned}

X_{ad} and X_{aq} are assumed to be determined by the generator’s magnetic saturation. The saturation characteristic is partitioned into three zones. The saturation characteristic is determined by open circuit.

7.3.1 A linear air gap region

The system p.u. gives this region a slope of one.

7.3.2 A nonlinear knee region

The knee region is model by exponential function V_{a} − V = A_{sat}e^{[B_{sat}(V −V_{L})]} The knee reactance X_{s} is model by

\begin{aligned}X_{s} & = & X_{u} \frac{V}{V_{a}} \\ & = & \frac{X_{u}}{1+ \frac{A_{sat}e^{[B_{sat}(V −V_{L})]}}{V}} \end{aligned}

X_{u} - unsaturated magnetising reactance V_{L} - knee region starting voltage V_{a} - air gap voltage corresponding to V

V_{a} − V is not zero at V_{L}. Hence, discontinuity is introduced at V_{L}.

7.3.3 A linear heavily saturated region

The saturated reactance X_{s} is model by \begin{aligned}X_{s} & = & X_{u} \frac{V}{V_{a}} \\ & = & X_{u} \frac{V}{R_{s}(V −V_{m})+V_{ma}} \end{aligned}

V_{ma} = V_{m} + A_{sat}e^{[B_{sat}(V −V_{L})]}

X_{u} - unsaturated magnetising reactance R_{s} - slope for heavily saturated region V_{m} - knee region ending voltage

7.4 Load Saturation

When under load, the saturation is approximated from the open circuit data.

\begin{aligned}Ψ_{a} & = & \sqrt{Ψ_{ad}^{2} + Ψ_{aq}^{2}}q \\ Ψ_{ad} & = & Ψ_{d} + X_{l}i_{d} \\ Ψ_{aq} & = & Ψ_{q} + X_{l}i_{q} \end{aligned}

7.5 Interfacing with the Network

In order to interface the generator mode to the rest of the network, we need to transform the actual voltages to synchronous rotating frame.

[ \begin{matrix}V_{D} \\ V_{Q} \end{matrix} ] = [ \begin{matrix}cos δ & − sin δ \\ sin δ & cos δ \end{matrix} ] [ \begin{matrix}V_{d} \\ V_{q} \end{matrix} ]
[ \begin{matrix}i_{D} \\ i_{Q} \end{matrix} ] = [ \begin{matrix}cos δ & − sin δ \\ sin δ & cos δ \end{matrix} ] [ \begin{matrix}i_{d} \\ i_{q} \end{matrix} ]

The transformation, T, is orthonormal. Hence, the inverse transformation T^{t}. i_{D}, i_{Q}, V_{D} and V_{Q} is determined by load-flow calculation. In directly, i_{d}, i_{q}, V_{d} and V_{q} is also determined by load-flow calculation. The steady state load angle is given by

tan(β + δ) = \frac{X_{q}P −r_{a}Q}{|V |^{2}+r_{a}P +X_{q}Q}

tan β = \frac{V_{D}}{V_{Q}} | V |^{2} = V_{D}^{2} + V_{Q}^{2} P = generated active power Q = generated reactive power

8 EPRI’s ETMSP Classical Synchronous Machine Model

Schematic of the classical generator model used by ETMSP is shown in figure 3

\begin{aligned}\frac{dω_{m}}{dt} & = & \frac{ω_{0}}{2H} (T_{m} − V_{D}i_{D} − V_{Q}i_{D}) \\ \frac{dδ}{dt} & = & ω_{m} \end{aligned}

8.1 Interfacing with the Network

\begin{aligned}V_{D} & = & E′ cos δ + X′_{d}I_{Q} \\ V_{Q} & = & E′ sin δ + X′_{d}I_{D} \end{aligned}

V_{D} and V_{Q} is determined by load-flow calculation. Indirectly, E′ is also determined by load-flow calculation and it is held constant over the cause of simulation.

Figure 3: EPRI detail synchronous machine model 9