Created by prashant iyengar last modified 2015-04-12

dcmotor1.py

      %%% DC motor first order model analysis
clear;

clear;

model=2;
if model==1
    J = 6e-7; % moment of inertia kg.cm2
    b = 0     ; % FRICTIONAL DAMPING
    Kt = 23.90e-3 ; %torque constant Nm/A
    Kb = 23.90e-3;  %back emf constant mV/rpm
    R = 10    ;%rotor resistance ohms
    L = 0.8e-3; %rotor inductance

    Tf=0.38e-3; %frrictional torque
    io=16e-03; % no load current

elseif model==0
    J = 6e-7; % moment of inertia kg.cm2
    b = 0     ; % FRICTIONAL DAMPING
    Kt = 12.30e-3 ; %torque constant Nm/A
    Kb = 12.30e-3;  %back emf constant mV/rpm
    R = 2.5    ;%rotor resistance ohms
    L = 0.22e-3; %rotor inductance

    Tf=0.38e-3; %frrictional torque
    io=31e-03; % no load current

elseif model==2
    J = 6e-7; % moment of inertia kg.cm2
    b = 0     ; % FRICTIONAL DAMPING
    Kt = 25.80e-3 ; %torque constant Nm/A
    Kb = 25.80e-3;  %back emf constant mV/rpm
    R = 15.2    ;%rotor resistance ohms
    L = 1e-3; %rotor inductance

    Tf=0.38e-3; %frrictional torque
    io=15e-03; % no load current

    elseif model==3
    J = 6e-7; % moment of inertia kg.cm2
    b = 0     ; % FRICTIONAL DAMPING
    Kt = 33.50e-3 ; %torque constant Nm/A
    Kb = 33.50e-3;  %back emf constant mV/rpm
    R = 32    ;%rotor resistance ohms
    L = 1.50e-3; %rotor inductance

    Tf=0.38e-3; %frrictional torque
    io=20e-03; % no load current

end
gain=12;

%% Second order DC motor model
s = tf('s');
P_motor = Kt/((J*s+b)*(L*s+R)+Kt*Kb);
disp('\n\nDC motor Step response parameters');
%step response
S = stepinfo(P_motor)
ltiview('step', P_motor, 0:0.1:5);
%[n,d]=tfdata(P_motor)
[a1,b1,c1]=zpkdata(P_motor);
dom=max(b1{1});
%p=poles(sys(P_motor));
disp(sprintf('Dominant Poles of DC motor  %0.03f',dom));

%% First order DC motor model
P_first = (Kt/R)/(J*s+b+(Kt*Kb/R));
disp('\n\nFirst order model DC motor Step response parameters');
%step response
S = stepinfo(P_first)
pole1=(b+(Kt*Kb/R))/J;
disp(sprintf('Poles of first DC motor model %0.03f  ',pole1));

ltiview('step', P_first, 0:0.1:5);

v=0:0.1:24;
angle=(v-R*io)*(Kt/(R))/((b+(Kt*Kb/R)));
figure(1)

xmarkers = [1,6,12,15,24]; % place markers at these x-values
ymarkers = (xmarkers-R*io)*(Kt/(R))/((b+(Kt*Kb/R)));


y1=ymarkers*60/(2*3.1412);
%plot(v,angle*60/(2*3.1412),'r')
plot(v,angle*60/(2*3.1412),'r',xmarkers,y1,'b*')
hold on

strValues = strtrim(cellstr(num2str([xmarkers(:) y1(:)],'(%d,%d)')));
text(xmarkers,y1,strValues,'VerticalAlignment','bottom');


grid on
xlabel('\bf voltage');
ylabel('\bf RPM');
title('\bf Voltage vs RPM');
hold off

    

dcmotor2.py

      clear;
model=1;
if model==1
    J = 6e-7; % moment of inertia kg.cm2
    b = 0     ; % FRICTIONAL DAMPING
    Kt = 23.90e-3 ; %torque constant Nm/A
    Kb = 23.90e-3;  %back emf constant mV/rpm
    R = 10    ;%rotor resistance ohms
    L = 0.8e-3; %rotor inductance

    Tf=0.38e-3; %frrictional torque
    io=16e-03; % no load current
    gain=12;
elseif model==0
    J = 6e-7; % moment of inertia kg.cm2
    b = 0     ; % FRICTIONAL DAMPING
    Kt = 12.30e-3 ; %torque constant Nm/A
    Kb = 12.30e-3;  %back emf constant mV/rpm
    R = 2.5    ;%rotor resistance ohms
    L = 0.22e-3; %rotor inductance

    Tf=0.38e-3; %frrictional torque
    io=31e-03; % no load current
    gain=6;
elseif model==2
    J = 6e-7; % moment of inertia kg.cm2
    b = 0     ; % FRICTIONAL DAMPING
    Kt = 25.80e-3 ; %torque constant Nm/A
    Kb = 25.80e-3;  %back emf constant mV/rpm
    R = 15.2    ;%rotor resistance ohms
    L = 1e-3; %rotor inductance

    Tf=0.38e-3; %frrictional torque
    io=15e-03; % no load current
    gain=15;
    elseif model==3
    J = 6e-7; % moment of inertia kg.cm2
    b = 0     ; % FRICTIONAL DAMPING
    Kt = 33.50e-3 ; %torque constant Nm/A
    Kb = 33.50e-3;  %back emf constant mV/rpm
    R = 32    ;%rotor resistance ohms
    L = 1.50e-3; %rotor inductance

    Tf=0.38e-3; %frrictional torque
    io=20e-03; % no load current
    gain=24;
end

%% Second order DC motor model
s = tf('s');
P_motor = Kt/((J*s+b)*(L*s+R)+Kt*Kb);
disp('\n\nDC motor Step response parameters');
S = stepinfo(P_motor)
%ltiview('step', P_motor, 0:0.1:5);
%[n,d]=tfdata(P_motor)
[a2,b1,c1]=zpkdata(P_motor);
dom=max(b1{1});
%p=poles(sys(P_motor));
disp(sprintf('Dominant Poles of DC motor  %0.03f',dom));

%% First order DC motor model
P_first = (Kt/R)/(J*s+b+(Kt*Kb/R));
%p=poles(sys(P_first));
%disp(sprintf('Poles of first DC motor model %0.03f  ',p(0)));
disp('\n\nFirst order model DC motor Step response parameters');
S = stepinfo(P_first)
pole1=(b+(Kt*Kb/R))/J;
disp(sprintf('Poles of first DC motor model %0.03f  ',pole1));

%ltiview('step', P_first, 0:0.1:5);

%% Voltage vs Speed
v=0:0.1:24;
angle=(v-R*io)*(Kt/(R))/((b+(Kt*Kb/R)));

figure(1);
xmarkers = [1,6,12,15,24]; % place markers at these x-values
ymarkers = (xmarkers-R*io)*(Kt/(R))/((b+(Kt*Kb/R)));


y1=ymarkers*60/(2*3.1412);
%plot(v,angle*60/(2*3.1412),'r')

plot(v,angle*60/(2*3.1412),'m',xmarkers,y1,'b*')
hold on
strValues = strtrim(cellstr(num2str([xmarkers(:) y1(:)],'(%d,%d)')));
text(xmarkers,y1,strValues,'VerticalAlignment','bottom');

grid on
xlabel('\bf voltage');
ylabel('\bf RPM');
title('\bf Voltage vs RPM');
hold off


%% torque speed characteristics
%Tf=0;
Tff=0;
Vs = gain-R*io;
t=0:0.05:30;
theta=((Vs*Kt/R)-(t/1000)-Tff)/(b+(Kb*Kt/R));
theta1=theta*60/(2*3.1412);
theta=theta1(theta1>0);
t=t(theta1>0);
figure(2);
power=theta.*((t/1000)+Tff)*(1000*2*3.1412/60);



%% Stall torque,current and no load speed computation

Tff=0;
stallTorque=(Vs*Kt/R)-Tff;
stallCurrent=(Vs/R);
OmegaNL = ((Vs*Kt/(R))-Tff)/( b + (Kt*Kt)/R);

fprintf('stall torque %0.3f\n',stallTorque);

fprintf('stall current %0.3f\n',stallCurrent);

fprintf('no load velocity %0.3f\n',OmegaNL*60/(2*3.1412));

hold on
Y=[0;stallTorque/2;stallTorque/2;0];
X=[0;0;OmegaNL/2;OmegaNL/2];

XX=X;%(index([1 1:end end]));
YY=Y;%[0 Y(index) 0];
      



Kt=(stallTorque*R/gain);
fprintf('computed value of Kt %0.3f\n',Kt);
b=(stallTorque/OmegaNL)-(Kt*Kt/R);
fprintf('computed value of b %0.3f\n',b);


%%%% motor efficiency computation

theta=theta*2*3.1412/60;
I = (Vs-(Kt*theta))/R;
J = (Kt*I) - (Kt*io)-Tff; 
Po = (J.*theta);
Pi = gain*I;
Eff = (Po./Pi)*100;
%figure(3);


[xx,yy]=max(Eff);
maxEff=xx;
mostEffVelocity=theta(yy);
disp(sprintf('Max Motor Efficiency : %f Percent ',100*maxEff));
disp(sprintf('Velocity at max efficiency %f rads/sec ',mostEffVelocity*60/(2*3.1412)));



%% plotting the motor characteristic curves
figure(2)
fill(XX*60/(2*3.1412),YY*1000,'b');
%


theta=theta*60/(2*3.1412);
[hAx,hLine1,hLine2] = plotyy(theta,t,theta,power);
Eff=Eff/2;
title ('\bf DC Motor Torque-Speed Characteristic'); 
xlabel('Speed')
ylabel(hAx(1),'Torque mNm') % left y-axis
ylabel(hAx(2),'Power mW') % right y-axis
h3=plot(theta,Eff,'r');
h4=plot(angle*60/(2*3.1412),v,'m',y1,xmarkers,'b*')

hold on
strValues = strtrim(cellstr(num2str([y1(:) xmarkers(:)],'(%d,%d)')));
text(y1,xmarkers,strValues,'VerticalAlignment','bottom');

legend([hLine1;hLine2;h3(1);h4(1)],'torque','power','efficiency','voltage');
%legend(h3(1),'efficiency')

%legend(h4(1),'voltage')

%axis([0 OmegaNL*60/(2*3.1412) 0 100]);
grid on;


%o = theta;
a11 = Kt*Vs/R-Kt*io;
b1 = ((Kt*Kt/R));
c1 = gain*Vs/R;
d1 = gain*Kt/R;
mostEffVelocity = ((b1*c1) - sqrt((b1*b1*c1*c1)-(a11*b1*c1*d1)))/(d1*b1);
disp(sprintf('Theorotical Velocity at max efficiency %f rads/sec ',mostEffVelocity*60/(2*3.1412)));


    

Comments (4)

HTTPS

SSH

You can clone a snippet to your computer for local editing. Learn more.