寒假期间预习李有堂编著、国防工业出版社出版的《机械系统动力学》,对其中第四章:单自由度系统的振动,第五章:两自由度系统的振动,的部分内容进行matlab编程仿真验证,并做简单研究。

1.1 有阻尼系统自由振动

2.1 谐波激励下的强迫振动2.2 阻尼对强迫振动的影响2.3 阻尼减振器

有阻尼系统自由振动在欠阻尼、临界阻尼和过阻尼三种情况下的位移曲线和相图如下图所示。

实现代码如下。

1%frictional oscillation simulation 2%lack friciton, critical friction and over friction 3%2019.02.19 4%Author: Jing-Xuan Yang 5%page: 73 6 7zeta = [3.4641, 34.641, 173.2051];%three conditions of friction 8tspan = linspace(0,2,400);%time 9%first condition10[t,x] = ode45('FreeOcillation',tspan,[0 1],[],zeta(1));%solve 2-ode11subplot(2,1,1);12plot(t,x(:,1),'-b');13hold on14subplot(2,1,2)15plot(x(:,1),x(:,2),'-b');16hold on17%second condition18[t,x] = ode45('FreeOcillation',tspan,[0 1],[],zeta(2));19subplot(2,1,1);20plot(t,x(:,1),'--g');21hold on22subplot(2,1,2)23plot(x(:,1),x(:,2),'--g');24hold on25%third condition26[t,x] = ode45('FreeOcillation',tspan,[0 1],[],zeta(3));27subplot(2,1,1);28plot(t,x(:,1),':r');29hold on30subplot(2,1,2)31plot(x(:,1),x(:,2),':r');32hold on33%draw titles and legends34subplot(2,1,1)35xlabel('Time(\tau)');36ylabel('Displacement x(\tau)');37title('Displacement as a function of \tau');38legend('\alpha = 0.1','\alpha = 1','\alpha = 5');39subplot(2,1,2);40xlabel('Displacement x(\tau)');41ylabel('Velocity v(\tau)');42title('Phase portrait');43legend('\alpha = 0.1','\alpha = 1','\alpha = 5');4445function xdot = FreeOcillation(t,x,dummy,zeta)4647xdot = [x(2);-2.0*zeta*x(2)-1200*x(1)];4849end1.2 谐波激励下的有阻尼强迫振动谐波激励下有阻尼强迫振动的幅频特性曲线和相频特性曲线如下图所示。

实现代码如下。

1%Phase Diagram for Different Damping Ratio 2%Amplitude-Frequency Characteristic Curve 3%2019/01/17 4%Author: Jing-Xuan Yang 5%page: 84-85 6 7x = 0:0.01:3; 8phi = ones(size(x)); 9H = ones(7,301);1011for i = 1:30112    ksi = 1;13    phi(i) = atan2(2*ksi*x(i),(1-(x(i))^2));14    H(1,i) = 1/sqrt((1-(x(i))^2)^2+(2*ksi*x(i))^2);15end16%Phase Diagram for Different Damping Ratio17figure(1)18plot(x,phi,'b.-');19hold on;2021for i = 1:30122    ksi = 0.5;23    phi(i) = atan2(2*ksi*x(i),(1-(x(i))^2));24    H(2,i) = 1/sqrt((1-(x(i))^2)^2+(2*ksi*x(i))^2);25end26plot(x,phi,'ko-');27hold on;2829for i = 1:30130    ksi = 0.25;31    phi(i) = atan2(2*ksi*x(i),(1-(x(i))^2));32    H(3,i) = 1/sqrt((1-(x(i))^2)^2+(2*ksi*x(i))^2);33end34plot(x,phi,'rx-');35hold on;3637for i = 1:30138    ksi = 0.15;39    phi(i) = atan2(2*ksi*x(i),(1-(x(i))^2));40    H(4,i) = 1/sqrt((1-(x(i))^2)^2+(2*ksi*x(i))^2);41end42plot(x,phi,'g+-');43hold on;4445for i = 1:30146    ksi = 0.10;47    phi(i) = atan2(2*ksi*x(i),(1-(x(i))^2));48    H(5,i) = 1/sqrt((1-(x(i))^2)^2+(2*ksi*x(i))^2);49end50plot(x,phi,'c*-');51hold on;5253for i = 1:30154    ksi = 0.05;55    phi(i) = atan2(2*ksi*x(i),(1-(x(i))^2));56    H(6,i) = 1/sqrt((1-(x(i))^2)^2+(2*ksi*x(i))^2);57end58plot(x,phi,'ys-');59hold on;6061for i = 1:30162    ksi = 0.00;63    phi(i) = atan2(2*ksi*x(i),(1-(x(i))^2));64    H(7,i) = 1/sqrt((1-(x(i))^2)^2+(2*ksi*x(i))^2);65end66plot(x,phi,'bp-');6768legend('\xi=1.00','\xi=0.50','\xi=0.25','\xi=0.15','\xi=0.10','\xi=0.05','\xi=0.00');69xlabel('\omega/\omega_n');70ylabel('\phi(\omega)');71title('Phase Diagram for Different Damping Ratio');7273%Amplitude-Frequency Characteristic Curve74figure(2)75plot(x,H(1,:),'b-');76hold on;77plot(x,H(2,:),'k-');78hold on;79plot(x,H(3,:),'r-');80hold on;81plot(x,H(4,:),'g-');82hold on;83plot(x,H(5,:),'c-');84hold on;85plot(x,H(6,:),'y-');86% hold on;87% plot(x,H(7,:),'b-');%goes to infinity88legend('\xi=1.00','\xi=0.50','\xi=0.25','\xi=0.15','\xi=0.10','\xi=0.05');89xlabel('\omega/\omega_n');90ylabel('|H(\omega)|');91title('Amplitude-Frequency Characteristic Curve');1.3 强迫振动的应用1.3.1 隔振原理主动隔振系数如下图所示。

测振原理如下图所示。

隔振原理以及测振原理的实现代码如下。

1%Active Vibration Isolation Coefficient 2%2019.01.19 3%Author: Jing-Xuan Yang 4%page: 108 5 6x = 0:0.01:5; 7nv = size(x);%size of x, vector 8n = nv(2);%number of x, scalar 9T = ones(7,n);10Beta = ones(7,n);1112for i = 1:n13    ksi = 1;14    T(1,i) = sqrt(1+(2*ksi*x(i))^2)/sqrt((1-(x(i))^2)^2+(2*ksi*x(i))^2);15    Beta(1,i) =x(i)*x(i)/sqrt((1-(x(i))^2)^2+(2*ksi*x(i))^2);16end1718for i = 1:n19    ksi = 0.5;20    T(2,i) = sqrt(1+(2*ksi*x(i))^2)/sqrt((1-(x(i))^2)^2+(2*ksi*x(i))^2);21    Beta(2,i) =x(i)*x(i)/sqrt((1-(x(i))^2)^2+(2*ksi*x(i))^2);22end2324for i = 1:n25    ksi = 0.25;26    T(3,i) = sqrt(1+(2*ksi*x(i))^2)/sqrt((1-(x(i))^2)^2+(2*ksi*x(i))^2);27    Beta(3,i) =x(i)*x(i)/sqrt((1-(x(i))^2)^2+(2*ksi*x(i))^2);28end2930for i = 1:n31    ksi = 0.15;32    T(4,i) = sqrt(1+(2*ksi*x(i))^2)/sqrt((1-(x(i))^2)^2+(2*ksi*x(i))^2);33    Beta(4,i) =x(i)*x(i)/sqrt((1-(x(i))^2)^2+(2*ksi*x(i))^2);34end3536for i = 1:n37    ksi = 0.10;38    T(5,i) = sqrt(1+(2*ksi*x(i))^2)/sqrt((1-(x(i))^2)^2+(2*ksi*x(i))^2);39    Beta(5,i) =x(i)*x(i)/sqrt((1-(x(i))^2)^2+(2*ksi*x(i))^2);40end4142for i = 1:n43    ksi = 0.05;44    T(6,i) = sqrt(1+(2*ksi*x(i))^2)/sqrt((1-(x(i))^2)^2+(2*ksi*x(i))^2);45    Beta(6,i) =x(i)*x(i)/sqrt((1-(x(i))^2)^2+(2*ksi*x(i))^2);46endfor i = 1:n50    ksi = 0.00;51    T(7,i) = sqrt(1+(2*ksi*x(i))^2)/sqrt((1-(x(i))^2)^2+(2*ksi*x(i))^2);52    Beta(7,i) =x(i)*x(i)/sqrt((1-(x(i))^2)^2+(2*ksi*x(i))^2);53end5455%Active Vibration Isolation Coefficient56figure(1)57plot(x,T(1,:),'b-');58hold on;59plot(x,T(2,:),'k-');60hold on;61plot(x,T(3,:),'r-');62hold on;63plot(x,T(4,:),'g-');64hold on;65plot(x,T(5,:),'c-');66hold on;67plot(x,T(6,:),'y-');68% hold on;69% plot(x,H(7,:),'b-');%goes to infinity70legend('\xi=1.00','\xi=0.50','\xi=0.25','\xi=0.15','\xi=0.10','\xi=0.05');71xlabel('\omega/\omega_n');72ylabel('|T(\omega)|');73title('Active Vibration Isolation Coefficient');7475%Vibration Measuring Device's Amplitude Diagram76figure(2)77plot(x,Beta(1,:),'b-');78hold on;79plot(x,Beta(2,:),'k-');80hold on;81plot(x,Beta(3,:),'r-');82hold on;83plot(x,Beta(4,:),'g-');84hold on;85plot(x,Beta(5,:),'c-');86hold on;87plot(x,Beta(6,:),'y-');88% hold on;89% plot(x,H(7,:),'b-');%goes to infinity90legend('\xi=1.00','\xi=0.50','\xi=0.25','\xi=0.15','\xi=0.10','\xi=0.05');91xlabel('\omega/\omega_n');92ylabel('|Beta(\omega)|');93title('Vibration Measuring Device''s Amplitude Diagram');941.3.3 惯性振动筛阿基米德螺线如下图所示。

实现代码如下。

1%Archimedes spiral 2%2019/01/19 3%Author: Jing-Xuan Yang 4 5t = 0:0.01:5; 6nv = size(t); 7n = nv(2); 8x = ones(1,n); 9y = ones(1,n);1011for i = 1:n12    x(i) = 3*t(i)*sin(5*t(i));13    y(i) = 3*t(i)*cos(5*t(i));14end1516plot(x,y,'-b');17xlabel('x');18ylabel('y');19title('Archimedes spiral');惯性振动筛的振幅与转速的关系如下图所示。

实现代码如下。

1%Inertial vibrating screen 2%2019/01/19 3%Author: Jing-Xuan Yang 4%page: 115 5 6omega = 0:0.01:5; 7nv = size(omega); 8n = nv(2); 9AM = ones(1,n);     %Amplitude10AMk = ones(1,n);    %Amplitude with different k11BM = ones(1,n);     %line12CM = ones(1,n);     %line13k1 = 45;14k2 = 5*k1;1516for i = 1:n17    AM(i) = abs(5*15*(omega(i))^2/(k1-(5+15)*(omega(i))^2));18    AMk(i) = abs(5*15*(omega(i))^2/(k2-(5+15)*(omega(i))^2));19    BM(i) = 15;20    CM(i) = 15/4;21end2223plot(omega,AM,'-b');24xlabel('\omega');25ylabel('|A_M|');26title('Inertial vibrating screen');27axis([0,5,0,100]);28text(-0.05,15,'r','Color','r');    %place text, absolute axis29text(-0.2,5,'|A''_M|','Color','c');30hold on;31plot(omega,AMk,'-g');32hold on;33plot(omega,BM,'-r');34hold on;35plot(omega,CM,'--c');2.1 谐波激励下的强迫振动两自由度系统谐波激励下强迫振动的幅频特性曲线如下图所示。

实现代码如下。

1%Two Degree of Freedom System Amplitude Diagram 2%2019/01/22 3%Author: Jing-Xuan Yang 4%page: 147 5 6omega = 0:0.01:5; 7omegan1 = 2; 8omegan2 = 4; 9nv = size(omega);10n = nv(2);11X1 = ones(1,n);12X2 = ones(1,n);13C = 10; %F/3k, a constant1415for i = 1:n16    X1(i) = C*(2-(omega(i)/omegan1)^2)/((1-(omega(i)/omegan1)^2)*(1-(omega(i)/omegan2)^2));17    X2(i) = C*1/((1-(omega(i)/omegan1)^2)*(1-(omega(i)/omegan2)^2));18end1920%Two Degree of Freedom System Amplitude Diagram21plot(omega,X1,'-b');22hold on;23plot(omega,X2,'--r');24legend('X_1(\omega)','X_2(\omega)');25xlabel('\omega');26ylabel('|X(\omega)|');27title('Two Degree of Freedom System Amplitude Diagram');2.2 阻尼对强迫振动的影响考虑阻尼时,两自由度系统谐波激励下强迫振动的幅频特性曲线如下图所示。

实现代码如下。

1%Two Degree of Freedom System Amplitude Diagram with Friction 2%2019/01/22 3%Author: Jing-Xuan Yang 4%Symbols: refer to <Dynamics of Mechanical System> Page 149 5 6m1 = 100; 7m2 = m1/20; 8mu = m2/m1; 9omega = 2:0.01:6;10omegan1 = 4;11omegan2 = 4;12alpha = omegan1/omegan2;13deltast = 1;1415lambda = omega/omegan1;16nlambda = size(lambda);17n = nlambda(2);1819Xi = [0 0.05 0.10 0.20 0.32 ];20nXi = size(Xi);21m = nXi(2);2223X1 = ones(m,n);2425for k = 1:m %friction2627    for i = 1:n %omega28        xi = Xi(k);29        numerator = ((lambda(i))^2-alpha^2)^2 + (2*xi*lambda(i))^2;30        denominator = (mu*(lambda(i))^2*alpha^2 - ((lambda(i))^2-1)*((lambda(i))^2-alpha^2))^2 ...31                       + (2*xi*lambda(i))^2*((lambda(i))^2-1+mu*lambda(i))^2;32        X1(k,i) = (1/deltast)*sqrt(numerator/denominator);33    end3435end3637%Two Degree of Freedom System Amplitude Diagram with Friction38plot(lambda,X1(1,:),'--b');39hold on;40plot(lambda,X1(2,:),'-r');41hold on;42plot(lambda,X1(3,:),'-y');43hold on;44plot(lambda,X1(4,:),'-c');45hold on;46plot(lambda,X1(5,:),'-g');47hold on;48plot(lambda,X1(6,:),'--k');49legend('\xi = 0','\xi = 0.05','\xi = 0.10','\xi = 0.20','\xi = 0.32','\xi = \infty');50xlabel('\lambda = \omega/\omega_{n1}');51ylabel('X_1/\delta_{st}');52title('Two Degree of Freedom System Amplitude Diagram with Friction');53axis([0.5 1.5 0 20]);54text(0.925,9,'S');55text(1.06,4.776,'T');5657%output abscissa for S and T, solving 4th polynomial function 58S = roots([1,mu/2, - (alpha^2*mu + 2*alpha^2 + 2)/2, - (alpha^2*mu)/2, alpha^2]);5960for j = 1:461    if S(j)>062        fprintf('Stationary points: lambda = %.4f\n',S(j));63    end64end在上面的幅频特性曲线中有两个不动点:S和T,求解此两个不动点的解析表达式的程序如下。

1%Solve two stationary points 2%2019/01/23 3%Author: Jing-Xuan Yang 4 5syms mu alpha lambda xi; 6 7numerator = (lambda)^2-alpha^2; 8denominator = (mu*(lambda)^2*alpha^2 - ((lambda)^2-1)*((lambda)^2-alpha^2))/ ... 9                       ((lambda)^2-1+mu*lambda);1011f = numerator == denominator ;1213solve(f, lambda)事实证明,用matlab进行符号运算求解表达式并不是很明智的选择。在这里,方程为四次代数方程。其求根公式异常复杂,故matlab只给出解是这个方程的根,而并未给出具体结果。对于解这个方程来说还不如手算,分类讨论即可。Matlab的可视化也不是很好,推荐使用maple软件做解各种方程(组)的符号运算,其公式的表达式可以与LaTeX无缝衔接。

阻尼减振器的类幅频特性曲线如下图所示。

课本上此图有一点错误,从上图可以看出这个图中也存在不动点,但是课本并未画出。本程序计算了该不动点的理论值,并与寻找交点算法得到的不动点值做了精确度的误差比较,实现代码如下。