1 Answers
Here is the code:
e=1.05;v=1.0;xdd=0.3;xl=0.4;pe0=1.0;h=5.0;tmax=4.0;dt=0.001;
f=60.0;d=0.5;tc=0.35;
ws=2*pi*f;
pmax=e*v/(xdd+0.5*xl);
pmech=pe0;
del=asin(pe0/pmax);
w=ws;
delta=[];time=[];n=1;
for t=0:dt:tmax
if t<=tc
xa=xdd+0.5*d*xl;xb=0.5*(1-d)*xl;xc=0.5*d*(1-d)*xl;
pmax=e*v*xc/(xa*xb+xb*xc+xc*xa);
else
pmax=e*v/(xdd+xl);
end
s_del=w-ws;
s_w=(pi*f/h)*(pmech-pmax*sin(del));
del=del+s_del*dt;
w=w+s_w*dt;
delta(n)=del*180/pi;
time(n)=t;
n=n+1;
end
plot(time,delta)