Chapter 8: Variable-Reluctance Machines and Stepping Motors

Example 8.1, Page number: 411

In [2]:
from __future__ import division
%matplotlib inline
from numpy import *
from math import *

#Variable declaration:
R=0.038                               #m
a=b=pi/3                            #rad
g=2.54*10**-4                       #m
D=0.13                              #m
N=100                               #turns in both poles
uo=4*pi*10**-7                      #permeability of free space()
i1=5                                #coil current (A)


#Calculation:
Lm=N**2*uo*a*R*D/(2*g)
#x=symbols('x')
subplot(2,1,1)
x=linspace(-180,-120,100)
L=-(Lm/60)*x-2*Lm
plot(x,L,'b')
#grid()

x=linspace(-60,0,100)
L=(Lm/60)*x+Lm
plot(x,L,'b')
grid()

x=linspace(0,60,100)
L=-(Lm/60)*x+Lm
plot(x,L,'b')
grid()


x=linspace(120,180,100)
L=(Lm/60)*x-2*Lm
plot(x,L)
annotate('Lm=0.128 H',xy=(-150,0.10))
annotate('Lmax',xy=(0,Lm+0.005))
ylabel('L11(theta)')
xlabel('theta')
grid()

#part(b)
subplot(2,1,2)
x1=linspace(-180,-120,100)
x2=linspace(-150,-90,100)
i1=5
i2=4
Tm1=(Lm/(2*pi/3))*i1**2
Tm2=(Lm/(2*pi/3))*i2**2
dll=np.ones(100)
plot(x1,-Tm1*np.array(dll),'g')
plot(x2,Tm2*np.array(dll),'b--')

x1=linspace(-60,0,100)
x2=linspace(-90,-30,100)
Tm1=(Lm/(2*pi/3))*i1**2
Tm2=(Lm/(2*pi/3))*i2**2
dll=np.ones(100)
plot(x1,Tm1*np.array(dll),'g')
plot(x2,-Tm2*np.array(dll),'b--')

x1=linspace(0,60,100)
x2=linspace(30,90,100)
Tm1=(Lm/(2*pi/3))*i1**2
Tm2=(Lm/(2*pi/3))*i2**2
dll=np.ones(100)
plot(x1,-Tm1*np.array(dll),'g')
plot(x2,Tm2*np.array(dll),'b--')

x1=linspace(120,180,100)
x2=linspace(90,150,100)
Tm1=(Lm/(2*pi/3))*i1**2
Tm2=(Lm/(2*pi/3))*i2**2
dll=np.ones(100)
plot(x1,Tm1*np.array(dll),'g')
plot(x2,-Tm2*np.array(dll),'b--')
grid()
ylim(-3,3)
annotate('___ i1=I1, i2=0', xy=(110,2.6))
annotate('---- i1=0, i2=I2', xy=(110,2.2))
ylabel('Torque [N.m]')
xlabel('thetam [degrees]')

#Results:
print "Lm =",Lm,"H"
print "(c)The peak torque =",round(Tm1,2),"N.m"
print "\t(i)  The net torque, (at thetam=0) =", 0, "N.m"
print "\t(ii) The net torque, (at thetam=45 deg.) =", 0, "N.m"
print "\t(iii)The net torque, (at thetam=75 deg) =", round(Tm1,2), "N.m"
show()

Populating the interactive namespace from numpy and matplotlib
Lm = 0.127968099059 H
(c)The peak torque = 1.53 N.m
	(i)  The net torque, (at thetam=0) = 0 N.m
	(ii) The net torque, (at thetam=45 deg.) = 0 N.m
	(iii)The net torque, (at thetam=75 deg) = 1.53 N.m

Example 8.3, Page number: 424

In [3]:
from __future__ import division
%matplotlib inline
from sympy import *
from math import *

#Variable declaration:
n1=4000                             #r/min
R=0.038                             #m
a=b=pi/3                            #rad
g=2.54*10**-4                       #m
D=0.13                              #m
N=100                               #turns in both poles
uo=4*pi*10**-7                      #permeability of free space(H/m)
Ll=0.005                            #H
Vo=100                              #phase voltage applied to phase 1.(V)


#Calculation:
wm=n1*pi/30
Lm=N**2*uo*a*R*D/(2*g)
thetam=symbols('thetam')
t=symbols('t')
#for part (a):
#for -60<=thetam<=0deg,
L11=Ll+(Lm/(pi/3))*(thetam+pi/3)
L111=diff(L11,thetam)
R1=L111*wm
#which is nuch greater than resistance R=1.5 ohm
thetam=-pi/3+wm*t
i1=Vo*t/(float(round(Ll,3))+float(Lm/(pi/3))*thetam+float(Lm/(pi/3))*pi/3)

#for part (b):
V2=-200                                     #applied voltage(V)
thetam2=symbols('thetam2')
L12=Ll+(Lm/(pi/3))*(pi/3-thetam2)
L112=diff(L12,thetam2)
to=2.5*10**-3                                     #ms
thetam2=float(-pi/3+wm*to)
i1=Vo*t/(float(round(Ll,3))+float(Lm/(pi/3))*thetam+float(Lm/(pi/3))*pi/3)
i2=(0.25-200*(t-to))/(0.005+51.1*(5*10**-3-t))


#Results:
print "i1 =",i1,"\t, (where round(16.2934044186179*pi,2) = 51.1 )"
print "\ni2 =",i2,"\n"


#Calculations & Results:
#for part (c):
from __future__ import division
from pylab import *

Lleak=0.005
Posintegral=0
integral=0
N1=500
tmax=3.75*10**-3
t=[0]*503
thet=[0]*503
Torque=[0]*503
deltat = tmax/N1
thetm=[0]*503
i=[0]*503
for n in range(1,N1+2,1):
    t[n-1]=tmax*(n-1)/N1
    thetm[n-1]=-(pi/3)+(400*pi/3)*t[n-1]
    if (thetm[n-1]<=0):
        i[n-1]=100*t[n-1]/(0.005+51.1*t[n-1])
        dld1d1theta = 0.122
        Torque[n-1]=0.5*i[n-1]**2*dld1d1theta
        Posintegral=Posintegral+Torque[n-1]*deltat
        integral=Posintegral
    else:
        i[n-1]=(0.25-200*(t[n-1]-2.5*10**-3))/(0.005+51.1*(5*10**-3-t[n-1]))
        dld11dtheta = -0.122
        Torque[n-1] = 0.5*i[n-1]**2*dld11dtheta
        integral = integral + Torque[n-1]*deltat

print "\nPositve torque integral =",Posintegral, "[N-m-sec]"
print "\nTorque integral=",integral,"[N-m-sec]\n"

plot(1000*np.array(t),i)
xlabel('time [msec]')
ylabel('Phase current [A]')
title('(a) phase-1 current profile')
grid()
show()
plot(1000*np.array(t),Torque)
xlabel('time [msec]')
ylabel('Torque [N-m]')
title('(b) torque profile')
grid()
show()

Populating the interactive namespace from numpy and matplotlib
i1 = 100*t/(51.1872396234976*t + 0.005) 	, (where round(16.2934044186179*pi,2) = 51.1 )

i2 = (-200*t + 0.75)/(-51.1*t + 0.2605) 


Positve torque integral = 0.000456384094483 [N-m-sec]

Torque integral= 0.000335463884625 [N-m-sec]


WARNING: pylab import has clobbered these variables: ['fmod', 'cosh', 'sinh', 'trunc', 'tan', 'gamma', 'degrees', 'radians', 'sin', 'expm1', 'ldexp', 'isnan', 'frexp', 'ceil', 'copysign', 'cos', 'tanh', 'fabs', 'sqrt', 'hypot', 'log', 'log10', 'pi', 'log1p', 'floor', 'modf', 'exp', 'isinf', 'e']
`%pylab --no-import-all` prevents importing * from pylab and numpy

Example 8.4, Page number: 433

In [4]:
from __future__ import division
%matplotlib inline

#Variavle declaration:
rpm=2500                        #rpm of motor


#Calculations & Results:
#For part (a):
theta=[0]*12
i=[0]*102
lambda1=[0]*102
for m in range(1,11,1):
    theta[m-1]=10*(m-1)
    for n in range(1,102,1):
        i[n-1]=30*(n-1)/100
        lambda1[n-1]=i[n-1]*(0.005+0.09*((90-theta[m-1])/90))*(8/(i[n-1]+8))

    
    plot(i,lambda1,'.')
    
    if m==1:
        hold(True)
        
xlabel('current [A]')
ylabel('Lambda [Wb]')
title('Family of lambda-i curves as theta_m varies from 0 to 90 degrees')    
annotate('theta_m=0 deg',xy=(6,0.03))
annotate('theta_m=0 deg',xy=(8,0.5))


#for part (b):
lambdamax=25*(0.005+0.09*(8/(25+8)))
AreaWnet=0
AreaWrec=0
deli=0.25
for n in range(1,102,1):
    i[n-1]=25*(n-1)/100
    AreaWnet=AreaWnet + deli*i[n-1]*(0.09)*(8/(i[n-1]+8))
    AreaWrec=AreaWrec + deli*(lambdamax-i[n-1]*(0.005+0.09*(8/(i[n-1]+8))))

Ratio=(AreaWnet+AreaWrec)/AreaWnet
print "part (b): Ratio =", round(Ratio,2)

#for part(b):
rps=rpm/60
T=1/rps
Pphase=2*AreaWnet/T
Ptot=2*Pphase
print "part (c): AreaWnet =", round(AreaWnet,2),"Joules"
print "Pphase =",round(Pphase),"W","\tPtot =",round(Ptot),"W\n"
plot(AreaWrec=0.7,AreaWnet=25)
grid()
show()

Populating the interactive namespace from numpy and matplotlib
part (b): Ratio = 1.55
part (c): AreaWnet = 9.91 Joules
Pphase = 825.0 W 	Ptot = 1651.0 W


WARNING: pylab import has clobbered these variables: ['power', 'random', 'fft', 'linalg', 'info']
`%pylab --no-import-all` prevents importing * from pylab and numpy