3-phase rectifier

In the diode rectifier given below, $v_a = 100\,\sqrt{2}\,\cos(100\,\pi\,t)$, $v_b = 100\,\sqrt{2}\,\cos(100\,\pi\,t - 90^{\circ})$, $v_c = 100\,\sqrt{2}\,\cos(100\,\pi\,t - 180^{\circ})$, $v_d = 100\,\sqrt{2}\,\cos(100\,\pi\,t - 270^{\circ})$. The magnitude of the current source is $10\,$A.
  1. Plot the voltage $v_o$.
  2. Determine the average value of $v_o$.
  3. What is the power delivered to $I_s$?
  4. Find the RMS values of the diode currents.
  5. Find the RMS value of the fundamental component of the ac line current.
  6. Find the input power factor.
  7. What is the THD of the input current?
  8. Which are the dominant harmonics in the ac line current?
In [1]:
from IPython.display import Image
Image(filename =r'rectifier_3ph_3_fig_1.png', width=360)
Out[1]:
No description has been provided for this image
In [2]:
# run this cell to view the circuit file.
%pycat rectifier_3ph_3_orig.in

We now replace the strings such as \$I0 with the values of our choice by running the python script given below. It takes an existing circuit file rectifier_3ph_3_orig.in and produces a new circuit file rectifier_3ph_3.in, after replacing \$I0 (etc) with values of our choice.

In [3]:
import gseim_calc as calc
import numpy as np

s_I0 = "10"

A_sin = 100.0*np.sqrt(2.0)

s_A_sin = ("%11.4E"%(A_sin)).strip()

l = [
  ('$I0', s_I0),
  ('$A_sin', s_A_sin),
]
calc.replace_strings_1("rectifier_3ph_3_orig.in", "rectifier_3ph_3.in", l)
print('rectifier_3ph_3.in is ready for execution')

rectifier_3ph_3.in is ready for execution
Execute the following cell to run GSEIM on rectifier_3ph_3.in.
In [4]:
import os
import dos_unix
# uncomment for windows:
#dos_unix.d2u("rectifier_3ph_3.in")
os.system('run_gseim rectifier_3ph_3.in')

get_lib_elements: filename gseim_aux/xbe.aux
get_lib_elements: filename gseim_aux/ebe.aux
Circuit: filename = rectifier_3ph_3.in
Circuit: n_xbeu_vr = 0
Circuit: n_ebeu_nd = 6
main: i_solve = 0
main: calling solve_trns
Transient simulation starts...
i=0
i=1000
i=2000
GSEIM: Program completed.
Out[4]:
0

The circuit file (rectifier_3ph_3.in) is created in the same directory as that used for launching Jupyter notebook. The last step (i.e., running GSEIM on rectifier_3ph_3.in) creates data files called rectifier_3ph_3_1.dat, etc. in the same directory. We can now use the python code below to compute/plot the various quantities of interest.

In [5]:
import numpy as np
import matplotlib.pyplot as plt 
import gseim_calc as calc
from setsize import set_size

f_hz = 50.0
T = 1.0/f_hz

slv = calc.slv("rectifier_3ph_3.in")

i_slv = 0
i_out = 0
filename = slv.l_filename_all[i_slv][i_out]
print('filename:', filename)
u1 = np.loadtxt(filename)
t1 = u1[:, 0]

col_v_o = slv.get_index(i_slv,i_out,"v_o")
col_ID1 = slv.get_index(i_slv,i_out,"ID1")
col_ID2 = slv.get_index(i_slv,i_out,"ID2")
col_ID3 = slv.get_index(i_slv,i_out,"ID3")
col_ID4 = slv.get_index(i_slv,i_out,"ID4")

i_out = 1
filename = slv.l_filename_all[i_slv][i_out]
print('filename:', filename)
u2 = np.loadtxt(filename)
t2 = u2[:, 0]

col_P_IS = slv.get_index(i_slv,i_out,"P_IS")
col_P_Vsa = slv.get_index(i_slv,i_out,"P_Vsa")
col_P_Vsb = slv.get_index(i_slv,i_out,"P_Vsb")
col_P_Vsc = slv.get_index(i_slv,i_out,"P_Vsc")
col_P_Vsd = slv.get_index(i_slv,i_out,"P_Vsd")

l_v_o   = calc.avg_rms_2(t1, u1[:,col_v_o  ], T, 2.0*T, 1.0e-5*T)
l_P_IS  = calc.avg_rms_2(t2, u2[:,col_P_IS ], T, 2.0*T, 1.0e-5*T)
l_P_Vsa = calc.avg_rms_2(t2, u2[:,col_P_Vsa], T, 2.0*T, 1.0e-5*T)
l_ID1   = calc.avg_rms_2(t1, u1[:,col_ID1  ], T, 2.0*T, 1.0e-5*T)

# get A_sin from the circuit file:
fin = open("rectifier_3ph_3.in", "r")
for line in fin:
    if 'name=Vsa' in line:
        for s in line.split():
            if s.startswith('a='):
                A_sin = float(s.split('=')[1])
fin.close()

print('average value of v_o:'          , "%11.4E"%( l_v_o  [1][0]))
print('rms value of v_o:'              , "%11.4E"%( l_v_o  [2][0]))
print('average power delivered to IS:' , "%11.4E"%(-l_P_IS [1][0]))
print('rms value of ID1:'              , "%11.4E"%( l_ID1  [2][0]))
print('average power delivered by Vsa:', "%11.4E"%( l_P_Vsa[1][0]))

Irms = l_ID1[2][0]
Vrms = A_sin/np.sqrt(2.0)
pf = l_P_Vsa[1][0]/(Vrms*Irms)

print('input power factor:', "%6.3f"%pf)

color1='red'
color2='goldenrod'
color3='blue'
color4='green'
color5='crimson'
color6='cornflowerblue'

fig, ax = plt.subplots(5, sharex=False)
plt.subplots_adjust(wspace=0, hspace=0.0)

set_size(5.5, 8, ax[0])

for i in range(5):
    ax[i].set_xlim(left=0.0, right=2.0*T*1e3)
    ax[i].grid(color='#CCCCCC', linestyle='solid', linewidth=0.5)

ax[0].set_ylabel(r'$v_o$',    fontsize=12)
ax[1].set_ylabel(r'$I_{D1}$', fontsize=12)
ax[2].set_ylabel(r'$I_{D2}$', fontsize=12)
ax[3].set_ylabel(r'$I_{D3}$', fontsize=12)
ax[4].set_ylabel(r'$I_{D4}$', fontsize=12)

ax[0].tick_params(labelbottom=False)
ax[1].tick_params(labelbottom=False)
ax[2].tick_params(labelbottom=False)
ax[3].tick_params(labelbottom=False)

ax[0].plot(t1*1e3, u1[:,col_v_o], color=color6, linewidth=1.0, label="$v_o$")

ax[1].plot(t1*1e3, u1[:,col_ID1], color=color1, linewidth=1.0, label="$I_{D1}$")
ax[2].plot(t1*1e3, u1[:,col_ID2], color=color2, linewidth=1.0, label="$I_{D2}$")
ax[3].plot(t1*1e3, u1[:,col_ID3], color=color3, linewidth=1.0, label="$I_{D3}$")
ax[4].plot(t1*1e3, u1[:,col_ID4], color=color4, linewidth=1.0, label="$I_{D4}$")

ax[4].set_xlabel('time (msec)', fontsize=12)

#plt.tight_layout()
plt.show()

filename: rectifier_3ph_3_1.dat
filename: rectifier_3ph_3_2.dat
average value of v_o:  1.2732E+02
rms value of v_o:  1.2793E+02
average power delivered to IS:  1.2732E+03
rms value of ID1:  4.9950E+00
average power delivered by Vsa:  3.1831E+02
input power factor:  0.637
No description has been provided for this image
In [6]:
import numpy as np
import matplotlib.pyplot as plt 
import gseim_calc as calc
from setsize import set_size

f_hz = 50.0
T = 1.0/f_hz

slv = calc.slv("rectifier_3ph_3.in")

i_slv = 0
i_out = 0
filename = slv.l_filename_all[i_slv][i_out]
print('filename:', filename)
u1 = np.loadtxt(filename)
t1 = u1[:, 0]

col_v_o = slv.get_index(i_slv,i_out,"v_o")
col_ID1 = slv.get_index(i_slv,i_out,"ID1")
col_ID2 = slv.get_index(i_slv,i_out,"ID2")
col_ID3 = slv.get_index(i_slv,i_out,"ID3")
col_ID4 = slv.get_index(i_slv,i_out,"ID4")

i_out = 1
filename = slv.l_filename_all[i_slv][i_out]
print('filename:', filename)
u2 = np.loadtxt(filename)
t2 = u2[:, 0]

col_P_IS = slv.get_index(i_slv,i_out,"P_IS")
col_P_Vsa = slv.get_index(i_slv,i_out,"P_Vsa")
col_P_Vsb = slv.get_index(i_slv,i_out,"P_Vsb")
col_P_Vsc = slv.get_index(i_slv,i_out,"P_Vsc")
col_P_Vsd = slv.get_index(i_slv,i_out,"P_Vsd")

# get A_sin from the circuit file:
fin = open("rectifier_3ph_3.in", "r")
for line in fin:
    if 'name=Vsa' in line:
        for s in line.split():
            if s.startswith('a='):
                A_sin = float(s.split('=')[1])
fin.close()

# compute Fourier coeffs:

t_start = T 
t_end = 2.0*T

n_fourier = 20

coeff_ID1, thd1_ID1, thd2_ID1 = calc.fourier_coeff_1A(t1, u1[:,col_ID1], 
    t_start, t_end, 1.0e-8, n_fourier)

coeff_v_o, thd_v_o = calc.fourier_coeff_1C(t1, u1[:,col_v_o], 
    t_start, t_end, 1.0e-8, n_fourier)

print("input current fundamental: RMS value: ", "%11.4E"%(coeff_ID1[1]/np.sqrt(2.0)))
print("input current THD (without dc term): ", "%11.4E"%thd1_ID1)
print("input current THD (with dc term): ", "%11.4E"%thd2_ID1)

x = np.linspace(0, n_fourier, n_fourier+1)

y_v_o = np.array(coeff_v_o)
y_ID1 = np.array(coeff_ID1)

fig, ax = plt.subplots(2, sharex=False)
plt.subplots_adjust(wspace=0, hspace=0.0)
grid_color='#CCCCCC'

set_size(6, 4, ax[0])

delta = 5.0
x_major_ticks = np.arange(0.0, (float(n_fourier+1)), delta)
x_minor_ticks = np.arange(0.0, (float(n_fourier+1)), 1.0)

for i in range(2):
    ax[i].set_xlim(left=-1.0, right=float(n_fourier))
    ax[i].set_xticks(x_major_ticks)
    ax[i].set_xticks(x_minor_ticks, minor=True)
    ax[i].grid(visible=True, which='major', axis='x', color=grid_color, linestyle='-', zorder=0)

ax[0].set_ylabel('$v_o$',fontsize=14)
ax[1].set_ylabel('$i_{D1}$',fontsize=14)

bars1 = ax[0].bar(x, y_v_o, width=0.3, color='red',   label="$v_o$", zorder=3)
bars2 = ax[1].bar(x, y_ID1  , width=0.3, color='green', label="$i_{D1}$" , zorder=3)

plt.tight_layout()
plt.show()

filename: rectifier_3ph_3_1.dat
filename: rectifier_3ph_3_2.dat
input current fundamental: RMS value:   3.1831E+00
input current THD (without dc term):   9.1958E-01
input current THD (with dc term):   1.2093E+00
No description has been provided for this image
In [7]:
import math
import numpy as np
import matplotlib.pyplot as plt 
import gseim_calc as calc
from setsize import set_size

f_hz = 50.0
T = 1.0/f_hz

slv = calc.slv("rectifier_3ph_3.in")

i_slv = 0
i_out = 0
filename = slv.l_filename_all[i_slv][i_out]
print('filename:', filename)
u1 = np.loadtxt(filename)
t1 = u1[:, 0]

col_v_o = slv.get_index(i_slv,i_out,"v_o")
col_ID1 = slv.get_index(i_slv,i_out,"ID1")
col_ID2 = slv.get_index(i_slv,i_out,"ID2")
col_ID3 = slv.get_index(i_slv,i_out,"ID3")
col_ID4 = slv.get_index(i_slv,i_out,"ID4")

# get A_sin from the circuit file:
fin = open("rectifier_3ph_3.in", "r")
for line in fin:
    if 'name=Vsa' in line:
        for s in line.split():
            if s.startswith('a='):
                A_sin = float(s.split('=')[1])
fin.close()

# compute Fourier coeffs:

t_start = T 
t_end = 2.0*T

n_fourier = 20

coeff_ID1, thd_ID1, coeff_a_ID1, coeff_b_ID1 = calc.fourier_coeff_2A(
  t1, u1[:,col_ID1], t_start, t_end, 1.0e-8, n_fourier)

t_ID1, y_ID1 = calc.construct_fourier_component(1, coeff_a_ID1, coeff_b_ID1, 0.0, T, 2, 200)

omg = f_hz*2.0*np.pi
vs_a = A_sin*np.sin(omg*t_ID1)

theta1 = math.atan2(-coeff_b_ID1[1], coeff_a_ID1[1])*180.0/math.pi

# if ID1 was written as k * sin (w*t + theta), what would theta be?
theta = theta1 + 90.0
print('angle of i_a w.r.t. V_a:', "%11.4E"%theta)

color1='green'
color2='crimson'
color3='cornflowerblue'

fig, ax = plt.subplots(2, sharex=False)
plt.subplots_adjust(wspace=0, hspace=0.0)

set_size(5.5, 5, ax[0])

for i in range(2):
    ax[i].set_xlim(left=0.0, right=2.0*T*1e3)
    ax[i].grid(color='#CCCCCC', linestyle='solid', linewidth=0.5)

ax[0].plot(t1*1e3, u1[:,col_ID1], color=color1, linewidth=1.0, label="$i_{D1}$")
ax[0].plot(t_ID1*1e3, y_ID1, color=color2, linewidth=1.0, linestyle='--', dashes=(5,3),
  label="$i_{D1}\,(1)$")

ax[1].plot(t_ID1*1e3, vs_a, color=color3, linewidth=1.0, label="$V_a$")

ax[1].set_xlabel('time (msec)', fontsize=12)

for i in range(2):
    ax[i].legend(loc = 'lower right',frameon = True, fontsize = 10, title = None,
      markerfirst = True, markerscale = 1.0, labelspacing = 0.5, columnspacing = 2.0,
      prop = {'size' : 12},)

plt.tight_layout()
plt.show()

filename: rectifier_3ph_3_1.dat
angle of i_a w.r.t. V_a:  1.4211E-14
No description has been provided for this image

This notebook was contributed by Prof. Nakul Narayanan K, Govt. Engineering College, Thrissur. He may be contacted at nakul@gectcr.ac.in.