《MATLAB计算方法和技巧62阻尼振动.docx》由会员分享，可在线阅读，更多相关《MATLAB计算方法和技巧62阻尼振动.docx（13页珍藏版）》请在三一办公上搜索。
1、MATLAB计算方法和技巧62阻尼振动弹簧振子的阻尼振动 问题一弹簧振子的质量为m，倔强系数为k。振子还受到与速度大小成正比、方向相反的阻力，比例系数为。当振子从静止开始运动时，初位移为A。物体的运动规律是什么？不同的阻尼下的运动曲线和速度曲线有什么特点？ 数学模型 根据牛顿运动定律，物体运动的微分方程为 m2dxdt22=-kx-gdxdt， (6.2.1) 取k/m = 0，/m = 2，0就是无阻尼时物体的固有角频率，是阻尼因子。物体的运动方程可表示为 dxdtrt22+2bdxdt+w0x=0。 (6.2.2) 2设微分方程的解为x = e，代入上式可得特征方程 r2 - 2r + 0
2、2 = 0。 (6.2.3) 特征方程的解为 r=-bb-w0， (6.2.4) 22设a=b-w0，可以是实数和零以及虚数，则r1 = - + ，r2 = - ，r1和r2可以是22实数或复数。微分方程的解为 x=C1e1+C2ertr2t=e-bt(C1eat+C2e-at)， (6.2.5) 其中C1和C2是由初始条件决定的常数。物体的速度为 v=dxdt=C1r1e1+C2r2ertr2t=e-btC1(-b+a)eat+C2(-b-a)e-at。 (6.2.6) 当t = 0时，x = A，v = 0，因此可得 A = C1 + C2，0 = C1(- + ) + C2(- - )，
3、 (6.2.7) 如果 0，即 0，解得两个常数分别为 a+ba-bC1=A，C2=A。 2a2a因此物体的位移为 x=A2ae-bt(a+b)eat+(a-b)e-at。 (6.2.8) 讨论当 0时，即 0，上式就是过阻尼的情况。 当 0时，即 0，不论用罗必塔法则还是用公式et 1 + t和e-t 1 - t，都可得 x=A(1+w0t)e-w0t。 (6.2.9) 这是临界阻尼的情况。 当 0时，设w=w0-b22，则 = i，i=-1为虚数单位，利用欧拉公式 1 e = cos + isin，e = cos - isin， 可得 cosq=12(eiqi-i+e-iq)，sinq=1
4、2i(eiq-e-iq)， 由(6.2.8)式可得 x=A2iwe-bt(iw+b)eiwt+(iw-b)e-iwt， 即 x=Ae-bt(cowst+或 x=其中j=arctan-bbw (6.2.10) swint， )w0wAe-btcosw(t+j，) (6.2.11) 。这就是欠阻尼的情况，物体可作准周期性运动，是其角频率，振幅w按指数规律衰减。物体作阻尼运动的周期为 T=2=2w， (6.2.12) 2w-b20可见：阻尼因子越大，周期越长。或者说：阻尼使振动变慢了。 当 = 0时，则 = i0，可得 x=A2eiw0t+e-iw0t=Acosw0t， 可见：在不计阻尼的情况下，物
5、体作简谐振动。 由于MATLAB能够进行复数运算，不论是过阻尼、临界阻尼，还是欠阻尼的情况，位移都可以用(6.2.8)式或(6.2.10)式计算。(6.2.8)式可用双曲函数表示 b-btx=Ae(coshat+sinhat)。 (6.2.13) a利用欧拉公式可以证明：coshi = cos，sinhi = isin，因此，(6.2.10)式和(6.2.13)式是完全相同的。 由(6.2.10)式可得物体运动的速度为 v=-Aw0w2e-btsinwt。 (6.2.14) 算法方法一：用解析解。取A为坐标的单位。取0的倒数为时间单位，则约化时间为t* = 0t；取约化阻尼因子为* = /0，
6、则约化阻尼角频率为 w=物体的坐标可表示为 *ww0=w0-bw022=1-b*2， x=exp(-bt)(coswt+*bw*sinwt)， (6.2.10*) *其中，x* = x/A。可见：阻尼因子* = /0决定了物体的位移曲线。取v0 = 0A为速度单位，物体的速度可表示为 2 v=*vv0=-1w*exp(-bt)sinwt。 (6.2.14*) 方法二：用微分方程的数值解。利用约化物理量，(6.2.2)式可化为 dxdt*2*2+2b*dxdt*+x=0。 (6.2.2*) *设x(1) = x，x(2) = dx/dt，可得 dx(1)dt*=x(2)，dx(2)dt*=-2b
7、x(2)-x(1)。 (6.2.2*) *由于 x(2)=dxdt*=1dxAw0dt=vAw0， 因此x(2)就是约化速度。初始条件为x(1) = 1，x(2) = 0。 方法三：用微分方程的符号解。将微分方程(6.2.2*)式可化为符号形式，利用dsolve指令可求符号解。 程序zqy6_2ode.m如下。 %阻尼运动的类型 clear %清除变量 b=input(请输入相对阻尼因子beta/w0:); %键盘输入阻尼因子 t=0:0.1:20; %固有角频率与时间的乘积w0t向量(约化时间向量) if b=1,b=b+eps;end %阻尼因子为1时则增加一小量 w=sqrt(1-b2)
8、; %准角频率向量 x=exp(-b*t).*(cos(w*t)+b/w*sin(w*t); %位移函数 v=-exp(-b*t).*sin(w*t)/w; %速度函数 figure %创建图形窗口 subplot(2,1,1) %选子图 plot(t,x) %画位移曲线 grid on %加网格 fs=16; %字体大小 xlabel(itomegarm_0itt,FontSize,fs)%标记横坐标 ylabel(itx/A,FontSize,fs) %标记纵坐标 title(质点阻尼运动的位移曲线,FontSize,fs)%标题 subplot(2,1,2) %选子图 plot(t,v)
9、 %画速度曲线 grid on %加网格 xlabel(itomegarm_0itt,FontSize,fs)%标记横坐标 ylabel(itv/omegarm_0itA,FontSize,fs)%标记纵坐标 title(质点阻尼运动的速度曲线,FontSize,fs)%标题 tm,X=ode45(zqy6_2fun,t,1;0,b);%计算位移和速度 hold on %保持图像 plot(t,X(:,2),r.) %画速度曲线 3 subplot(2,1,1) %选子图 hold on %保持图像 plot(t,X(:,1),r.) %画位移曲线 s=D2x+,num2str(2*b),*D
10、x+x; %微分方程字符串 sx=dsolve(s,x(0)=1,Dx(0)=0); %微分方程的符号解 sv=diff(sx); %求速度的符号解 x=subs(sx,t,t); %位移 v=subs(sv,t,t); %速度 plot(t,x,ko) %画位移曲线 legend(解析解,数值解,符号解) %图例 txt=itbeta/omegarm_0:,num2str(b);%参数文本 text(0,0,txt,FontSize,fs) %显示文本 subplot(2,1,2) %选子图 plot(t,v,ko) %画速度曲线 legend(解析解,数值解,符号解) %图例 程序要调用函数zqy6_2fun.m。 %阻尼运动的二阶微分方程的函数 function f=fun(t,x,flag,b) f=x(2);-2*b*x(2)-x(1); %速度和加速度 说明用解析解画了曲线之后，可取图片，按任一键继续。用微分方程的数值解和符号解所画的点和圈与公式法的曲线完全重合，说明三种方法都是正确的。 P6.2a图 P6.2b图 P6.2c图 作业：用微分方程的数值解计算和绘制单摆阻尼振动的曲线。 湖南大学物电院 周群益 4
