Short circuit current calculations in Python. Full article at: https://firsttimeprogrammer.blogspot.it/2016/11/short-circuit-currents-calculation-on.html

short_circuit_current.py

import matplotlib.pyplot as plt
from scipy.integrate import quad
import numpy as np
plt.style.use("ggplot")

Vn = 132000                             # Nominal voltage [V]
E = Vn * np.sqrt(2) / np.sqrt(3)        # Peak phase voltage [V]
l = 80                                  # Line length [km]
L = 0.00048 * l                         # Line inductance [H]
R = 0.0754 * l                          # Line resistance [Ohm]
F = 50                                  # Frequency [Hz]
w = 2*np.pi*F                           # Frequency [rad / s]
z = np.sqrt(R**2 + (w*L)**2)            # Impedance (magnitude) [Ohm]
tau = L/R                               # Time constant [s]
alpha_phi = -np.pi/2                    # Alpha - phi [rad]. Use worst case scenario.
N = 10000                               # Number of points used
t = np.linspace(0, 8*tau, N)            # Time [s]

# Unidirectional component of short circuit current (a)
def unidirectional(t, alpha_phi=alpha_phi, tau = tau, E=E, z=z):
    y = - (np.sqrt(2)*E/z)*np.exp(-t/tau) * np.sin(alpha_phi)
    return y

# Sinusoidal component of short circuit current (b)
def sinusoidal(t, alpha_phi=alpha_phi, E=E, z=z):
    y = (np.sqrt(2)*E/z)*np.sin(w*t+alpha_phi)
    return y

# Short circuit current = a + b
def short_cc_current(t):
    y = sinusoidal(t) + unidirectional(t)
    return y

# Short circuit current squared
def short_cc_current_squared(t):
    y = short_cc_current(t)**2
    return y

# Current RMS
Icc_rms = np.sqrt( (1/np.max(t)) * quad(short_cc_current_squared, 0, np.max(t))[0] )
# Short circuit current
Icc = sinusoidal(t) + unidirectional(t)

# Print info
print("Tau:", round(tau*1000, 5), "ms")
print("Alpha - phi:", round(alpha_phi/np.pi*180, 3),"degree")
print("Frequency:",F ,"Hz")
print("Line reactance:", round(w*L/l, 3),"Ohm/km")
print("Line resistance:", round(R/l, 3),"Ohm/km")
print("Line length:",l,"km")
print("Series impedance:", round(z, 3),"Ohm")
print("Impedance characteristic angle:", round(np.arctan((w*L)/R)*180/np.pi, 3), "degrees")
print("Nominal voltage rms:", round(Vn/1000, 3), "kV")
print("Phase voltage rms:", round(E/(np.sqrt(2)*1000), 3), "kV")
print("Peak positive Icc current:", round(np.max(Icc)/1000, 3), "kA")
print("Peak negative Icc current:", round(np.min(Icc)/1000, 3), "kA")
print("Icc rms", round(Icc_rms/1000, 3), "kA")

# Plot data
plt.figure()
plt.plot(t, Icc, color='blue', label = "Short circuit current")
plt.plot(t, [np.max(Icc)]*t.shape[0], color="red", label="max")
plt.plot(t, [np.min(Icc)]*t.shape[0], color="red", label="min")
plt.plot(t, unidirectional(t), color="green", label = "unidirectional comp")
plt.xlim([0, np.max(t)])
plt.ylim([np.min(Icc)-1000, np.max(Icc) + 8500])
plt.title("Icc")
plt.legend()
plt.xlabel("Time [t]")
plt.ylabel("Short circuit current [A]")

plt.show()

################################################################################
#                           OUTPUT
################################################################################

#Tau: 6.36605 ms
#Alpha - phi: -90.0 degree
#Frequency: 50 Hz
#Line reactance: 0.151 Ohm/km
#Line resistance: 0.075 Ohm/km
#Line length: 80 km
#Series impedance: 13.488 Ohm
#Impedance characteristic angle: 63.434 degrees
#Nominal voltage rms: 132.0 kV
#Phase voltage rms: 76.21 kV
#Peak positive Icc current: 13.714 kA
#Peak negative Icc current: -11.28 kA
#Icc rms 8.158 kA