%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)])