单自由度系统的响应图，相位图，根轨迹图绘制
function free_vibration_with_vious_damping(m,k,c,x0,dx0,t)%equation : M*ddx+C*dx+k*x=0
%again asuming : xt=c*exp(s*t),
%and then inserting this function into equation
%obtain the characteristic equation : m*s^2+c*s+k=0(the roots is s1 and s2
%these roots give two solution: xt1=C1*exp(s1*t) and xt2=C2*exp(s2*t)
%thus the general solution is : x(t)=C1*exp(s1*t)+C2*exp(s2*t)
% where C1and C2 is determined by the initial conditions of the system
%howerer S12=(-c+/-sqrt(C^2-4mk))/2M; (2)
%the radical in the equation(2)becomes zero as Critial Damping Cc=2*sqrt(k*m)
%damping radio e=C/Cc;
%then wei can rewrite the equation of (2) as :
% S12=(-e +/- sqrt(e^2-1))*wn=-e*wn +/- wd (wd=sqrt(1-e^2)*wn)
%now wei must discuss three situation
%Frist one named underddamped system (0<e<1)
% S1=(-e+i*sqrt(1-e^2))*wn and S2=(-e - i*sqrt(1-e^2))*wn
%at last X(t)=expt(-e*wn*t)*{ C1*cos(wd*t)+C2*sin(wd*t)}
% where C1=x0;C2=(dx0+e*wn*x0)/wd
%the second one e=1
% S1=S2=-wn
%because of the repeated rotos,the solution of equation of(1) is
%X(t)=(C1+C2*t)*exp(-wn*t)
% where C1=x0 ;C2=dx0+wn*x0
%the thrid situation overdamped system (e>1)
% S1=(-e+sqrt(e^2-1))*wnS2=(-e-sqrt(e^2-1))*wn
%X(t)=C1*exp(S1*t)+C1*exp(S2*t)
%where C1=(x0*wn(e+sqrt(e^2-1))+dx0)/(2*wn*sqrt(e^2-1))
% C2= (-x0*wn(e-sqrt(e^2-1))-dx0)/(2*wn*sqrt(e^2-1))
clc
clear all
%M=450;K=26519.2;C=1000; x0=0.539;dx0=1;t=0:0.001:10; %checking!!
M=4;K=4;C=20;x0=0.539;dx0=1;t=0:0.001:10;
=size(C);
%if nargin<6 t=0:0.01:20;end
if m==n==1
wn=sqrt(K/M); e=C/(2*sqrt(K*M)); %damping radio
if e<1&e>0
wd=sqrt(1-e^2)*wn;
C1=x0;C2=(dx0+e*wn*x0)/wd;
Xt=exp(-e*wn.*t).*( C1*cos(wd.*t)+C2*sin(wd.*t));
elseif e==1
C1=x0 ;C2=dx0+wn*x0;
Xt=(C1*ones(size(t))+C2.*t).*exp(-wn.*t);
elseif e>1
S1=(-e+sqrt(e^2-1))*wn ; S2=(-e-sqrt(e^2-1))*wn;
C1=(x0*wn*(e+sqrt(e^2-1))+dx0)/(2*wn*sqrt(e^2-1)) ;
C2= (-x0*wn*(e-sqrt(e^2-1))-dx0)/(2*wn*sqrt(e^2-1));
Xt=C1.*exp(S1.*t)+C1.*exp(S2.*t);
end
figure(1)
plot(t,Xt,'r'); grid on; xlabel('time');ylabel('diaplacement');title('free vibration with vious damping');
end
%plot the Locus of S1 and S2 of this system
wn=sqrt(K/M); Cc=2*sqrt(K*M);
e1=; e2=;
S1=(-e1+i*sqrt(1-e1.^2)).*wn;S2=(-e1 - i*sqrt(1-e1.^2)).*wn;
S3=(-e2+sqrt(e2.^2-1)).*wn ; S4=(-e2-sqrt(e2.^2-1)).*wn;
S1r=real(S1);S1i=imag(S1);S2r=real(S2);S2i=imag(S2);
S3r=real(S3);S3i=imag(S3);S4r=real(S4);S4i=imag(S4);
figure(2)
plot(S1r,S1i,'r');grid on; hold on;
plot(S2r,S2i,'b');title('locus of this systen');hold on
plot(S3r,S3i,'*'); plot(S4r,S4i,'-');xlabel('real axis');
%plot Phane plane of damped system 使用R-K方法
x1=x0;x2=dx0;T=0;h=0.001;n=1;x1old=x1;x2old=x2;
while T<=10*pi/wn;
K11=x2old;
K12=-(C*x2old+K*x1old)/M;
T=T+h/2; x1=x1old+ K11*h/2;x2=x2old+0.5*h*K12;
K21=x2;
K22=-(C*x2+K*x1)/M;
T=T+h/2;x1=x1old-h*K11+2*h*K21;x2=x2old-h*K12+2*h*K22;
K31=x2;
K32=-(C*x2+K*x1)/M;
x1old=x1old+(h/6)*(K11+4*K21+K31);
x2old=x2old+(h/6)*(K12+4*K22+K32);
x11(1,n)=x1old;
x21(1,n)=x2old;
n=n+1;
end
figure(3)
plot(x11,x21,'r');
grid on
xlabel('displacement')
ylabel('speed')
title('phase plane')
end
本帖最后由 ChaChing 于 2011-10-23 16:27 编辑
原创 or 转贴 !? 个人专业有限, 看似像分享!? 是吗?
若是, 若能稍作简易说明下相关输入/输出/算法, 甚至举例, 或许更适合学习! 个人浅见
页:
[1]