import math
#initialisation of variables
WR2 = 5.82
Smach = 1333
n = 1800
#Calculations
ft_lb = 746 / 550.0
w = 2 * math.pi * n / 60
H = ft_lb * WR2 * w**2 / (2 * 32.2 * Smach)
#Results
print('The inertia constant in MJ/MVA is = %.2f' %H)
print('Converting H to a 100-MVA system base,units in MJ/MVA = %.2f' %(H * Smach / 100))
import math
#initialisation of variables
P1 = 500.0
pf1 = 0.85
V1 = 20.0
n1 = 3600
P2 = 1333.0
pf2 = 0.9
V2 = 22.0
n2 = 1800.0
Pbase = 100
H1 = 4.8
H2 = 3.27
#Calculations
KE = H1 * P1 + H2 * P2
#Results
print('The total kinetic energy of rotation of the two machines in MJ is = %.2f' %KE)
print('The inertia constant for the equivalent machine on 100-MVA base in MJ/MVA is = %.2f' %(KE/Pbase))
import math
#initialisation of variables
Pm =1
Vt = 1.0
V_ib = 1.0
X1_g =0.2
X1_t = 0.1
X1_l1 = 0.4
X1_l2 = 0.4
#Calculations
X = X1_t + X1_l1 /2
a = math.asin(Pm*X/(Vt*V_ib))*180/math.pi
Vt1 = Vt*complex(math.cos(a * math.pi / 180),math.sin(a * math.pi / 180))
I = (Vt1 - V_ib) / complex(X)
E1 = Vt1 + complex(X1_g * I)
X1 = X1_g + X1_t + X1_l1 /2
Pmax = abs(E1) * V_ib / X1
#Results
print('The terminal voltage after considering the angle is,in per unit= {0:.2f}+{1:.2f}i'.format(Vt1.real, Vt1.imag))
print('The output current from the generator in per unit is= {0:.2f}+{1:.2f}i'.format(I.real, I.imag))
print('The transient internal voltage in per unit= {0:.2f}+{1:.2f}i'.format(E1.real, E1.imag))
print('Power angle equation is')
print("\n Pe = %.2f * sin(delta) \n where delta is the machine rotor angle wrt to the infinite bus" %Pmax)
import math
#initialisation of variables
H = 5.0
Pm =1.0
Vt = 1.0
V_ib = 1.0
X1_g =0.2
X1_t = 0.1
X1_l1 = 0.4
X1_l2 = 0.4
X = X1_t + X1_l1 /2
a = math.asin(Pm * X / (Vt * V_ib)) * 180 / math.pi
Vt1 = Vt*complex(math.cos(a * math.pi / 180),math.sin(a * math.pi / 180))
I = (Vt1 - V_ib) / complex(X)
E1 = Vt1 + complex(X1_g * I)
y10 = complex(3.33)
y32 = complex(2.5)
y30 = complex(5)
y20 = complex(5)
Ybus = [[0,0,0],[0,0,0],[0,0,0]]
#Calculations
Ybus[0][0] = -y10
Ybus[0][1] =0
Ybus[0][2] = y10
Ybus[1][0] = Ybus[0][1]
Ybus[1][1] = -(y32 + y30)
Ybus[1][2] = y32
Ybus[2][0] = Ybus[0][2]
Ybus[2][1] = Ybus[1][2]
Ybus[2][2] = -(y10 + y32 + y30)
print('Ybus formed by inspection is')
print(Ybus)
Ybus_new = [[0,0],[0,0]]
for c in range(0, 2):
for d in range(0, 2):
Ybus_new[c][d] = Ybus[c][d] - ((Ybus[c][2]*Ybus[2][d]) / Ybus[2][2])
print('Ybus formed after elimination of Bus 3')
print(Ybus_new)
Pmax = abs(E1) * V_ib * abs(Ybus_new[0][1])
delta = 28.44
Pa = Pm - Pmax * math.sin(delta * math.pi / 180)
b = 180 * Pa / H
#Results
print('The power abgle equation is')
print(" Pe = %.3f * sin(delta) \n where delta is the machine rotor angle wrt to the infinite bus" %Pmax)
print('The swing equation is')
print("(%.2f/180f) * d(delta)^2/dt^2 = %.2f - %.2fsin(delta) \n" %(H,Pm,Pmax))
print(" Intial Accelerating power is %.3f per unit \n" %(Pa))
print(" Initial acceleration is %.2f*f \n where f is the system frequency" %(b))
import math
#initialisation of variables
H = 5
Pm =1
Vt = 1
V_ib = 1
X1_g =0.2
X1_t = 0.1
X1_l1 = 0.4
X1_l2 = 0.4
#Calculations
X = X1_t + X1_l1 /2
a = math.asin(Pm * X / (Vt * V_ib)) * 180 / math.pi
Vt1 = Vt * complex(math.cos(a * math.pi / 180),math.sin(a * math.pi / 180))
I = (Vt1 - V_ib) / complex(X)
E1 = Vt1 + complex(X1_g * I)
y12 = 1 / complex((X1_g + X1_t + X1_l1))
Y12 = -y12
Pe = abs(E1) * V_ib * abs(Y12)
#Results
print('The power abgle equation is')
print(" Pe = %.3f * sin(delta) \n where delta is the machine rotor angle wrt to the infinite bus" %Pe)
print('The swing equation is')
print(" (%.2f/180f) * d(delta)^2/dt^2 = %.2f - %.2fsin(delta) \n" %(H,Pm,Pe))
import math
#initialisation of variables
delta = 28.44
H = 5
ws = 377
Pm =1
Vt = 1.0
V_ib = 1
X1_g =0.2
X1_t = 0.1
X1_l1 = 0.4
X1_l2 = 0.4
#Calculations
X = X1_t + X1_l1 /2
a = math.asin(Pm * X / (Vt * V_ib)) * 180 / math.pi
Vt1 = Vt * complex(math.cos(a * math.pi / 180),math.sin(a * math.pi / 180))
I = (Vt1 - V_ib) / complex(X)
E1 = Vt1 + complex(X1_g * I)
X1 = X1_g + X1_t + X1_l1 /2
Pmax = abs(E1) * V_ib / X1
Sp = Pmax * math.cos(delta * math.pi / 180)
wn = math.sqrt(ws * Sp / (2 * H))
fn = wn / (2 * math.pi)
T = 1 / fn
#Results
print("The angular frequency of oscillation is %.3f elec rad/s " %wn)
print("The corresponding frquency of oscillation is %.2f Hz " %fn)
print("The period of oscillation is %.3f s" %T)
import math
#initialisation of variables
delta = 28.44
H = 5
ws = 377
Pm =1
Vt = 1
V_ib = 1
X1_g =0.2
X1_t = 0.1
X1_l1 = 0.4
X1_l2 = 0.4
#Calculations
X = X1_t + X1_l1 /2
a = math.asin(Pm * X / (Vt * V_ib)) * 180 / math.pi
Vt1 = Vt * complex(math.cos(a * math.pi / 180),math.sin(a * math.pi / 180))
I = (Vt1 - V_ib) / complex(X)
E1 = Vt1 + complex(X1_g * I)
X1 = X1_g + X1_t + X1_l1 /2
Pmax = abs(E1) * V_ib / X1
delta_rad = delta * math.pi / 180
delta_cr = math.acos((math.pi - 2 * delta_rad) * math.sin(delta_rad) - math.cos(delta_rad))
t_cr = math.sqrt(4 * H * (delta_cr - delta_rad) / (ws * Pm))
#Results
print("Critical clearing angle = %.3f elec rad " %(delta_cr))
print("Critical clearing angle for the system = %.3f s" %(t_cr))