G = 50.0 #Rating of machine(MVA)
f = 50.0 #Frequency of turbo generator(Hz)
V = 11.0 #Voltage rating of machine(kV)
H = 9.0 #Cycle corresponding to 180 ms
P_0 = 40.0 #Pre-fault output power(MW)
delta_0 = 20.0 #Rotor angle at instant of fault(degree)
P_0_close = 0 #Output power at instant of reclosing(MW)
P_a = P_0 - P_0_close #Net accelerating power(MW)
delta_sqr = P_a*180*f/(G*H) #double derivative(elect.degrees/sec^2)
from scipy.integrate import quad
def integrand1(t): #Integrates the double derivative to 800*t
return delta_sqr
a, err = quad(integrand1, 0, 180*10**-3) #Rotor velocity(electrical degrees/sec)
def integrand2(t): #Integrates the double derivative to 400*t^2
return delta_sqr*t
b, err = quad(integrand2, 0, 180*10**-3)
delta = delta_0 + b #Rotor angle(electrical degrees)
print('Rotor angle at the instant of reclosure = %.2f° electrical degrees' %delta)
print('Rotor velocity at the instant of reclosure = %.1f electrical degrees/sec' %a)
V = 1.0 #Infinite bus voltage(p.u)
E = 1.0 #e.m.f of finite generator behind transient reactance(p.u)
X_T = 0.8 #Transfer reactance(p.u)
P_i = 0.5 #Input power(p.u)
P_i_d = 0.8 #p.u
P_0 = 0.5 #Output power(p.u)
P = 0.5 #Power(p.u)
import math
P_m = E*V/X_T #Amplitude of power angle curve(p.u)
delta_0 = math.asin(P_i/P_m) #Radians
delta = math.asin(P_i_d/P_m) #Radians
delta_m = math.pi-delta #Radians
A_acc = P_i_d*(delta-delta_0)-P_m*(math.cos(delta_0)-math.cos(delta)) #Possible area of accceleration
A_dec = P_m*(math.cos(delta)-math.cos(delta_m))-P_i_d*(delta_m-delta) #Possible area of deceleration
if(A_acc < A_dec):
print('System is stable')
stability = A_dec/A_acc
print('Margin of stability = %.2f' %stability)
else:
print('System is not stable')
x = 0.25 #Transient reactance(p.u)
E = 1.0 #e.m.f of finite generator behind transient reactance(p.u)
x_T = 0.1 #Reactance of transformer(p.u)
x_L = 0.4 #Reactance of one line(p.u)
P_i = 0.25 #Pre-fault power(p.u)
import math
X_T = x+x_T+(x_L/2) #Transfer reactance at pre-fault state(p.u)
P_m = E**2/X_T #Amplitude of power angle curve at pre-fault state(p.u)
X_T1 = 1.45 #Transfer reactance b/w finite generator & infinite bus at faulty state(p.u).Refer texbook problem for figure
P_m1 = E**2/X_T1 #Amplitude of power angle curve at faulty state(p.u)
r1 = X_T/X_T1
delta_0 = math.asin(P_i/P_m) #Radians
delta_1 = math.asin(P_i/(r1*P_m)) #Radians
delta_m = math.pi - delta_1 #Radians
from scipy.integrate import quad
def integrand1(delta):
return r1*P_m*math.sin(delta)
a, err = quad(integrand1, delta_0, delta_1)
A_acc = P_i*(delta_1-delta_0) - a
def integrand2(delta):
return r1*P_m*math.sin(delta)
b, err = quad(integrand2, delta_1, delta_m)
A_dec = b - P_i*(delta_m-delta_1)
limit = 0.5648 #Obtained by iterations.Refer textbook.Here assigned directly.
if(A_acc < A_dec):
print('System is Stable')
stability = A_dec/A_acc
print('Margin of stability = %.2f' %stability)
else:
print('System is not stable')
print('Transient stability limit = %.4f p.u' %limit)
print('\nNOTE : ERROR : angle delta_0 = 7.9° = 0.13788 radian not 0.014 radian as in textbook')
x = 0.25 #Transient reactance(p.u)
E = 1.0 #e.m.f of finite generator behind transient reactance(p.u)
x_T = 0.1 #Reactance of transformer(p.u)
x_L = 0.4 #Reactance of one line(p.u)
P_i = 0.7 #Pre-fault power(p.u)
import math
X_T = x+x_T+(x_L/2) #Transfer reactance at pre-fault state(p.u)
P_m = E**2/X_T #Amplitude of power angle curve at pre-fault state(p.u)
X_T1 = 1.45 #Transfer reactance b/w finite generator & infinite bus at faulty state(p.u).Refer texbook problem for figure
P_m1 = E**2/X_T1 #Amplitude of power angle curve at faulty state(p.u)
r1 = X_T/X_T1
X_T2 = x+x_T+x_L #Transfer reactance for post fault state(p.u)
r2 = X_T/X_T2
P_m2 = r2*P_m
delta_0 = math.asin(P_i/P_m) #Radians
delta_1 = math.asin(P_i/(r2*P_m)) #Radians
delta_m = math.pi - delta_1 #Radians
delta_c = 0.7 #Specified value(radians)
from scipy.integrate import quad
def integrand1(delta):
return r1*P_m*math.sin(delta)
a, err = quad(integrand1, delta_0, delta_c)
A_acc = P_i*(delta_c-delta_0) - a
def integrand2(delta):
return r2*P_m*math.sin(delta)
b, err = quad(integrand2, delta_c, delta_m)
A_dec = b - P_i*(delta_m-delta_c)
cos_delta_cr = ((delta_m-delta_0)*math.sin(delta_0)-r1*math.cos(delta_0)+r2*math.cos(delta_m))/(r2-r1)
delta_cr = math.acos(cos_delta_cr)*180/math.pi
if(A_acc < A_dec):
print('System is Stable')
stability = A_dec/A_acc
print('Margin of stability , K = %.2f' %stability)
else:
print('System is not stable')
print('Critical clearing angle for a certain pre-fault power = %.2f°' %delta_cr)
print('Critical clearing time will be known from circuit-breaker specifications')
P_i = 0.75 #Pre-fault power(p.u)
f = 50.0 #Frequency(Hz)
H = 6.0 #Value of H for finite machine(sec)
x_G = 0.2 #Reactance of machine(p.u)
x_T = 0.1 #Reactance of transformer(p.u)
x_L = 0.4 #Reactance of line(p.u)
V = 1.0 #Voltage of infinite bus(p.u)
E = 1.0 #e.m.f of finite generator behind transient reactance(p.u)
import math
X_T = x_G+x_T+(x_L) #Transfer reactance at pre-fault state(p.u)
P_m = E**2/X_T #Amplitude of power angle curve at pre-fault state(p.u)
delta_0 = math.asin(P_i/P_m) #Radians
delta_0a = delta_0*180/math.pi
delta_cr = math.acos((math.pi-2*delta_0)*math.sin(delta_0)-math.cos(delta_0))
delta_cra = delta_cr*180/math.pi
t_cr = ((delta_cra-delta_0a)*2*H/(180*f*P_i))**0.5
print('Critical clearing angle for circuit breaker at bus 1 = %.2f°' %delta_cra)
print('Time for circuit breaker at bus 1 ,t_cr = %.3f sec' %t_cr)