Vh=30*0.44704;B=2200;K2=10440;K1=5410;t=[0, 20];M=345;m=20;
opts=odeset('MaxStep',0.1);
[T,X]=ode23s(@(t,x)p5fun(t,x,Vh,B,K1,K2),[0 20],[0,0,0,0],opts);
y=X(:,1);z=X(:,2);vy=X(:,3);vz=X(:,4);
figure(1)
subplot(2,2,1);plot(T,y)
subplot(2,2,2);plot(T,z)
subplot(2,2,3);plot(T,vy)
subplot(2,2,4);plot(T,vz)
dotvy=(K2/M)*(z-y)+(B/M)*(vz-vy);
figure(2)
plot(T,dotvy);
ymax=max(abs(y));
Arms=sqrt((1/20)*trapz(T,dotvy.^2));
fprintf('The max vertical displacement is %g metres. \n',ymax)
fprintf('The root mean square of acceleration is %g m/s^2.\n',Arms)
Vht=[10:5:80]*0.44704;
for k=1:length(Vht)
[T1,X1]=ode23s(@(t,x)p5fun(t,x,Vht(k),B,K1,K2),[0 20],[0,0,0,0],opts);
[ymax1(k),Arms1(k)]=p5_3fun(Vh,B,K1,K2);
end
figure(3);plot(Vht,ymax1)
figure(4);plot(Vht,Arms1)
B1o=[1000:200:4000];
K2o=(5000:1000:15000);
for p=1:length(K2o)
[T1,X1]=ode23s(@(t,x)p5fun(t,x,Vh,B,K1,K2o(p)),[0 20],[0,0,0,0],opts);
[ymax2(p),Arms2(p)]=p5_3fun(Vh,B,K1,K2o(p));
end
figure(5); plot(K2o,ymax2)
figure(6);plot(K2o,Arms2)
for q=1:length(B1o)
[T1,X1]=ode23s(@(t,x)p5fun(t,x,Vh,B,K1,K2),[0 20],[0,0,0,0],opts);
[ymax3(q),Arms3(q)]=p5_3fun(Vh,B1o(q),K1,K2);
end
figure(7);plot(B1o,ymax3)
figure(8);plot(B1o,Arms3)
B1o=[1000:200:4000];
K2o=(5000:1000:15000);
for p=1:lgenth(K2o)
for q=1:length(B1o)
[T1,X1]=ode23s(@(t,x)p5fun(t,x,Vh,B1o(q),K1,K2(p),[0 20],[0,0,0,0],opts));
[ymax4(q,p),Arms4(q,p)]=p5_3fun(Vh,B1o(q),K1,K2o(p));
end
end
figure(9);mesh(K2o,B1o,ymax4)
figure(10);mesh(K2o,B1o,Arms4)
function u=roadprofile(L)
if (L>2) && (L<=2.5)
u=-0.2*(L-2);
elseif (L>2.5) && (L<=3.5)
u=-0.1;
elseif (L>3.5) && (L<=4)
u=0.2*(L+4);
elseif (L>4) && (L<=100)
u=0;
elseif (L>100) && (L<=104)
u=-0.15*sin((pi/4)*L-(28*pi));
else
u=0;
end
end
function dotx=p5fun(t,x,Vh,B,K1,K2)
L=Vh*t; u=roadprofile(L); M=345;m=20;
dotx=[x(3);x(4);
(K2/M)*(x(2)-x(1))+(B/M)*(x(4)-x(3));
(K1/m)*(u-x(2))-(K2/m)*(x(2)-x(1))-(B/m)*(x(4)-x(3))];
end
function [ymax,Arms]=p5_3fun(Vht,B,K1,K2)
global X1 T1
M=345;
dotvy1=(K2/M).*(X1(:,2)-X1(:,1))+(B/M).*(X1(:,4)-X1(:,3));
ymax=max(abs(X1(:,1)));
Arms=sqrt((1/20)*trapz(T1,dotvy1.^2));
end
0 Comments
Sign in to comment.