%Motor2.m -A Ruina, Feb 16, 2010
%This is NOT adequately commented
disp('************************************')
disp(['START Calculation at ' datestr(now)])
clear all %added on 2/21/10
syms omega V I T k c R x D
%The basic motor equations
eq1 = T- k*I + c*omega;
eq2 = V- k*omega - I*R;
%Find I,V from T and omega
S=solve(eq1,eq2, I, V);
%Write power in and out in terms of T and omega
Pin = S.I*S.V; %S.I and S.V are in terms of T and omega
Pout = T*omega;
eff1 = simplify(Pout/Pin);
eff2 = subs(eff1, T, x*omega); % that is, x is defined to be T/omega
eff3=simplify(eff2)
%Use first derivative to find the local max
bestx = solve(diff(eff3, x)) %eff replaced by eff3 on 2/21/10
bestx = bestx(1) % just keep the positive root, added on 2/21/10
% Substitute back in to get best efficiency
besteff1 = simple(subs(eff3, x, bestx));
% Try to tidy up and simplify
besteff2=simple(subs(besteff1, k, D*sqrt(R*c))); % defines D = k/sqrt(R*c)
pretty(besteff2)
%Specifications for Faulhaver 2657 - 012 CR
k= 0.0173; %in consistant SI units
R=0.71; %ohms. actually, the spec is wrong and R = 1.3
c=3e-6; %in consistant SI units
D = k/sqrt(R*c) %The key motor parameter
besteffnumerical = subs(besteff2)
bestxnumerical = subs(bestx) % Remember, this is the ratio T/omega
%Motor use specifications, given by design team, say
Pout = 30; % in watts
omega_out = 5; % in rad/s
% We know the bestx = T/omega. We know Pout = T*omega. So
omegamotor = sqrt(Pout/bestxnumerical)
G = omegamotor/omega_out
Tmotornum = bestxnumerical * omegamotor
Vbat = subs(S.V); Vbat = subs(Vbat, [T omega],[Tmotornum omegamotor])
disp(['END Calculation at ' datestr(now)])