物理学阻尼振动中阻尼系数的研究 大学毕业设计.docx

1、物理学阻尼振动中阻尼系数的研究 大学毕业设计阻尼振动中阻尼系数的研究摘 要：在我们的力学教材和其他的有关力学知识的书中，对于阻尼振动中阻尼系数的研究已屡见不鲜，并且有好多学者和老师都对它进行了分析和研究，并取得了不错的理论结果。但在我们课本中只是利用了传统的试验方法来测定阻尼振动中的阻尼系数，在本文中，我们实验测得的数据进行拟合来计算空气阻尼和磁阻尼的阻尼系数，得出了令人满意的结果，与传统实验手段比较起来更为准确。关键词：阻尼系数；阻尼振动；空气阻尼；磁阻尼；曲线拟合A study on damping coefficient of vibration dampingABSTRACT: In

2、our mechanical materials and other books relating to the mechanical knowledge, for damping vibration damping constants have been frequent, and there are many scholars and teachers of his were analyzed and studied and made good theoretical results. But in our textbooks, just use the traditional test

3、methods to measure the damping in vibration damping, in this paper, we use curve fitting with the experimental data combined with the specific procedures in Matlab to calculate the air damping and magnetic damping, the damping coefficient, obtained satisfactory results, and experimental means a lot

4、more convenient.KEYWORDS: Damping coefficient; damping vibration; air damping; magnetic damping; curve fitting引 言机械振动是力学中的常见现象，实际生活中的振动并不是像理想条件下的简谐振动，而是带有阻尼的阻尼振动。最近研究阻尼对振动的影响也成为人们比较关注的问题。在实验上传统的测定方法是通过利用半衰期或利用阻尼振动曲线求出阻尼振子的阻尼系数。在本文中我们根据实验数据利用Matlab做出理论曲线并同时画出速度峰值的拟合曲线，再将二者进行比较估计的方法求得阻尼系数。在本文中，我们不仅利用这

5、种方法讨论只在空气阻力下的阻尼系数并且画出拟合曲线，而且我们还讨论了在磁阻尼作用下的阻尼系数。文中不仅讨论了阻尼系数在两种情况下与速度一次方的关系，还讨论了与速度二次方的关系。测量过程简单易行且得到的实验数据可靠, 具有一定的推广价值。在下面我们分别对仅存在空气阻尼及存在磁阻尼的情况进行讨论。1 数据的曲线拟合工程技术实践中, 特别是在实验测试中, 往往只能得到两个相关变量的一系列离散值, 它们间的函数关系就只得用列表法或图示法表示。但是，有时因为某种原因希望把这种关系转换成解析表达式。另外,有些变量间的关系虽然可以用解析法表示，由于数学解析式过于繁杂不便使用，也希望用一个简单的解析表达式来近

6、似地代替它，凡此种种，利用简单解析式近似地代替列表法、图示法或复杂解析式表示的函数一类问题，即寻找一个满足的简单函数的问题，都可称之为函数逼近。数据拟合法就是函数逼近的重要方法之一。科学实验或数据统计中，在得不到两组相关数据之间的精确解析表达式时，就希望找到关系的近似解析式，越简单越好。在这种情况下，我们提出了函数逼近方法-数据拟合法，它不要求构造的近似函数全部通过样本点，而是“很好的逼近”它们。这种近似函数反映了已知数据组间存在着某种关系的一般趋势，是“拟合”这些数据得出的函数曲线，不同于“插值”得出的曲线. 寻找“很好的逼近”的函数有多种方法，我们常用的是最小二乘法。在MatLab中，有专

7、用的拟合指令，指令是polyfit，使用格式为：p=polyfit(x,y,m)，可以方便地给出求算结果，我们便利用实验得到的振动数据结果利用拟合方法，对数据点进行拟合，得到拟合曲线。 2 模 型对于存在阻尼力的情况下，阻尼振动1方程为： (2.1)其中,、分别为系统的阻尼系数和固有频率. 方程(2.1)的解为： (2.2)其中，、分别为弹簧振子的振幅、初始相位角，而，为阻尼振动周期 (2.3)由方程(2.2)可得其速度与时间的关系式为： (2.4)这为我们的实验测量2提供了理论依据，根据公式所有公式都与有关，因此只需要记录时间，周期即可。3 无磁阻尼时空气阻尼的阻尼系数的研究实验所采用的气垫

8、导轨型号是L-QG-Y-1500/5.8，电子计数器型号是MUJ-5C。在气垫导轨上测量弹簧振子振动时的阻尼系数，进而求得时间常数及品质因数等物理量。假定振动为线性振动，且气垫导轨存在的空气阻尼很小。测量时所用滑块质量为179.70g，振动的振幅为30cm。由实验测得弹簧振子振动连续10个周期的数据，见表3-1。表3-1振动周期测量数据n12345678910T(s)1.26901.26911.26911.26911.26901.26891.26901.26881.26871.2687利用Matlab计算10个周期测量值的算术平均值、实验平均值的标准误差。计算程序如下：clearhn=1.26

9、90,1.2691,1.2691,1.2691,1.2690,1.2689,1.2690,1.2688,1.2687,1.2687; %填入10个周期测量值nn=size(hn);N=nn(2);hp=0;for n=1:N; hp=hp+hn(n)/N;endhpfor n=1:N; vhn(n)=abs(hn(n)-hp)2;endsigma2=0;for n=1:N; sigma2=sigma2+vhn(n)/(N-1);endsigman=sqrt(sigma2);s= sigman/sqrt(N);hp, s %endhp代表算术平均值, s为平均值的实验标准误差。输出结果：hp =

10、 1.2689s= 4.9889e-005即：. 在实验时取初始位移，初始速度为，并通过放置在平衡位置的电子计数器记录下连续10个周期中，一个周期从正方向和反方向通过两次平衡位置时光电门的速度。 表3-2滑块通过平衡位置时的速度0.251.252.253.254.255.256.257.258.259.25145.98142.45139.86137.36134.95132.28129.87127.71125.63123.300.751.752.753.754.755.756.757.758.759.75-143.88-141.04-138.50-136.05-133.69-131.06-128

11、.86-126.74-124.53-122.25程序设计：clearT=1.2689;te=0.25*T,1.25*T,2.25*T,3.25*T,4.25*T,5.25*T,6.25*T,7.25*T,8.25*T,9.25*T; %输入时间ve=145.98,142.45,139.86,137.36,134.95,132.28,129.87,127.71,125.63,123.30; %输入速度pv=polyfit(te,ve,2); %生成拟合多项式(向量形式),阶次为3pvt1=poly2str(pv,t); %生成拟合多项式(字符形式)pvt1输出结果：pvt1 =0.032136

12、t2 - 2.3349 t + 146.4178clearT=1.2689;te=0.75*T,1.75*T,2.75*T,3.75*T,4.75*T,5.75*T,6.75*T,7.75*T,8.75*T,9.75*T; %输入时间ve=-143.88,-141.04,-138.50,-136.05,-133.69,-131.06,-128.86,-126.74,-124.53,-122.25; %输入速度pv=polyfit(te,ve,2); %生成拟合多项式(向量形式),阶次为2pvt2=poly2str(pv,t); %生成拟合多项式(字符形式)pvt2输出结果：pvt2 =-0.0

13、21055 t2 + 2.1588 t - 145.8269，分别是沿正方向运动时，通过平衡位置的光电门的时间和速度。 ，是沿负方向运动时，通过平衡位置的光电门的时间和速度。 3.1阻尼系数由阻尼振动方程的解(2.4)式编制Matlab程序计算理论上时间与速度曲线, 同时将两拟合曲线绘制在同一幅图中。程序设计如下：clearA=30;phi=-pi;deta=0.017;T=1.2689;omegaf=2*pi/T;t=0:0.001:14;vn1=0.032136*t.2- 2.3349*t+146.4178;vn2=-0.021055*t.2+2.1588*t-145.8269;v=-de

14、ta*A*exp(-deta*t).*cos(omegaf*t+phi)-A*omegaf*exp(-deta*t).*sin(omegaf*t+phi);plot(t,vn1,-.k,t,vn2,:k,t,v,-k,LineWidth, 2)xlabel (fontsize14rmt(T), Color,k)axis(0 10*T -200 200)set(gca,XTick,0:T:10*T) %给定x轴的刻度set(gca,XTickLabel,0,1,2,3,4,5,6,7,8,9,10) %给定x轴每个刻度的标示值ylabel (fontsize14rmv(cm/s), Color,

15、k)legend(拟合曲线1, 拟合曲线2, 理论曲线)grid on%end该程序绘制出弹簧振子时间与速度曲线, 见图2-1与2-2。其中实线、点线是由实验上得到向正方向、负方向通过平衡位置光电门的时间与速度数据拟合的曲线, 实线是由理论方程(2.3)所得到的曲线3。即为所要求的阻尼系数，在用Matlab画图时的值是在不断地调试的，当将解中的调整为时, 理论上曲线通过平衡位置时的速度与两拟合曲线基本吻合, 即两按拟合曲线正是速度曲线的包络线。所以得到此时阻尼系数为： (4)可以看出通过对的多次取值估计，重复多次计算画图，我们得出了较为满意阻尼系数，使得理论线与拟合曲线基本吻合，如图2-1所示

16、，结果证明上述理论与实验结合的方法是可行的。与此同时，我们还可以得到时间常数和品质因数： 时间常数为 (5) 品质因数为 (6)但是我们从图2-2看出拟合曲线并没有很好的与理论实线相吻合，放大后观察到还有不相交的部分，也就是不完全包络，原因在于我们绘制的拟合曲线是一次直线，所以结果并不是令人很满意。 ；无论对于怎么调整，对于放大的图像来看，总不可能全部包络速度的峰值，只能是尽可能的平衡，这说明仅仅调整是不够的，还要考虑阻尼与速度的二次方的关系。3.2 考虑阻尼与速度二次方有关上面所讨论的是空气阻尼与速度一次方的拟合关系，在某些点上并不能很好的拟合，下面利用MatLab来模拟与速度二次方5的拟合

17、曲线。在M文档中建立如下函数：function Q=funxz(t,x)T=1.2689;omegaf=2*pi/T;d=0.0145;p=0.00043;Q=x(2);-omegaf2*x(1)-2*d*x(2)-2*p*x(2).2;程序如下：cleart,x=ode45(funxz,0:0.01:14,-30 0);T=1.2689;t1=0:0.001:14;vn1=0.032136*t1.2- 2.3349*t1+146.4178;vn2=-0.021055*t1.2+2.1588*t1-145.8269;plot(t1,vn1,-k,t1,vn2,-k,t,x(:,2),-k,Li

18、neWidth, 2)xlabel (fontsize14rmt(T), Color,k)axis(0 10*T -200 200)set(gca,XTick,0:T:10*T) %给定x轴的刻度set(gca,XTickLabel,0,1,2,3,4,5,6,7,8,9,10) %给定x轴每个刻度的标示值ylabel (fontsize14rmv(cm/s), Color,k)legend(拟合曲线1, 拟合曲线2, 数值计算曲线)grid on%end在M文档中所建立的函数恰好是阻尼系数与速度二次方的函数关系式，我们又对和p进行了多次调试，结果如图2-3与图2-4所示，通过对阻尼与速度二次

19、方关系的讨论，所得结果分别与图2-1和图2-2比较，拟合曲线与速度峰值更好的吻合了，并且得到了更为精确的，并且得出了速度二次方因子p=0.00043。；。4 研究磁阻尼情况下的阻尼系数在滑块两侧对称地粘贴上永久小磁铁块，此时由于磁阻尼4力的作用，使振动的振幅衰减的更快。在此时测量的阻尼系数，仍假定振动为线性振动。测量时所用滑块质量为191.06g，振动的振幅为30.00cm。由实验测得弹簧振子连续振动10个周期的数据，见表4-1。表4-1振动周期测量数据n12345678910T(s)1.30001.30001.29961.29981.29961.29991.29991.29981.29971

20、.2998利用MatLab计算10个周期测量值的算术平均值、实验平均值的标准误差。计算程序：clearhn=1.3000,1.3000,1.2996,1.2998,1.2996,1.2999,1.2999,1.2998,1.2997,1.2998; %填入10个周期测量值nn=size(hn);N=nn(2);hp=0;for n=1:N; hp=hp+hn(n)/N;endfor n=1:N; vhn(n)=abs(hn(n)-hp);endfor n=1:N; vhn2(n)=(hn(n)-hp)2;endsigma2=0;for n=1:N; sigma2=sigma2+vhn2(n)/

21、(N-1);endsigman=sqrt(sigma2);s= sigman/sqrt(N);hp, s %hp代表算术平均值, s为平均值的实验标准误差. %endhp =1.2998s=4.5826e-005即：。在实验时取初始位移，初始速度为。我们启记录下连续10个周期中，一个周期从正方向和反方向通过两次平衡位置时光电门的速度。表4-2振动周期测量数据0.251.252.253.254.255.256.257.258.259.25140.65135.32131.06126.74122.85119.05114.94111.23107.76104.490.751.752.753.754.75

22、5.756.757.758.759.75-137.74-133.33-128.86-124.84-120.92-117.23-113.25-109.41-106.04-102.77，是向正方向运动时，通过平衡位置的光电门的时间和速度。根据上述数据编程拟合出，的关系式。程序设计：clearT=1.2998;te=0.25*T,1.25*T,2.25*T,3.25*T,4.25*T,5.25*T,6.25*T,7.25*T,8.25*T,9.25*T; %输入时间ve=140.65,135.32,131.06,126.74,122.85,119.05,114.94,111.23,107.76,10

23、4.49; %输入速度pv=polyfit(te,ve,2); %生成拟合多项式(向量形式),阶次为3pvt1=poly2str(pv,t); %生成拟合多项式(字符形式)pvt1输出结果：；由上述结果可得与之间的关系为：、是向负方向运动时，通过平衡位置的光电门的时间和速度。 程序设计：clearT=1.2998;te=0.75*T,1.75*T,2.75*T,3.75*T,4.75*T,5.75*T,6.75*T,7.75*T,8.75*T,9.75*T; %输入时间ve=-137.74,-133.33,-128.86,-124.84,-120.92,-117.23,-113.25,-109

24、.41,-106.04,-102.77; %输入速度pv=polyfit(te,ve,2); %生成拟合多项式(向量形式),阶次为3pvt2=poly2str(pv,t); %生成拟合多项式(字符形式)pvt2输出结果：。由上述结果可得与之间的关系为： 4.1 阻尼系数由阻尼振动方程(3)式的解编制Matlab程序计算理论上给出的时间与速度曲线，同时将两拟合曲线绘制在同一幅图中。程序中的vn1即，vn2即。程序设计：clearA=30;phi=-pi;deta=0.028;T=1.2998;omegaf=2*pi/T;t=0:0.001:14;vn1=0.049482*t.2-3.6731*t

25、+141.5108;vn2=-0.037352*t.2+3.5008*t-141.0681;v=-deta*A*exp(-deta*t).*cos(omegaf*t+phi)-A*omegaf*exp(-deta*t).*sin(omegaf*t+phi);plot(t,vn1, -.k,t,vn2,:k,t,v,-k,LineWidth, 2)xlabel (fontsize14rmt(T), Color,k)axis(0 10*T -200 200)set(gca,XTick,0:T:10*T) %给定x轴的刻度set(gca,XTickLabel,0,1,2,3,4,5,6,7,8,9,

26、10) %给定x轴每个刻度的标示值ylabel (fontsize14rmv(cm/s), Color,k)legend(拟合曲线1, 拟合曲线2, 理论曲线)grid on%end该程序绘制出弹簧振子时间与速度曲线，见图3-1。其中点划线、点线是由实验上得到向正方向、负方向通过平衡位置光电门的时间与速度数据拟合的曲线，实线是理论方程所得到的曲线。当将解中的调整为时，理论上曲线通过平衡位置时的速度与两拟合曲线基本吻合。由此得到阻尼系数为 (4) 时间常数为 (5) 品质因数为 (6)； 如图3-2.、图3-1所示，我们得到了基本吻合的拟合曲线，但放大后观察，不论如何调整，还是没有得到满意的效果

27、，那么就要考虑与速度二次方的关系。 4.2 考虑速度二次方因子上面所讨论的是磁阻尼与速度一次方的拟合关系，在某些点上并不能很好的拟合，下面利用MatLab来模拟与速度二次方5的拟合曲线。在M文档中建立如下函数：function Q=funxz(t,x)T=1.2998;omegaf=2*pi/T;d=0.025;b=0.00068;Q=x(2);-omegaf2*x(1)-2*d*x(2)-2*b*x(2).2;程序如下：cleart,x=ode45(funxz,0:0.01:14,-30 0); %t,x=ode45(funerx,0:0.01:2*pi,4 0);T=1.2998;t1=0

28、:0.01:14;vn1=0.049482*t1.2-3.6731*t1+141.5108;vn2=-0.037352*t1.2+3.5008*t1-141.0681;plot(t1,vn1,-k,t1,vn2,-k,t,x(:,2),-k,LineWidth, 2)xlabel (fontsize14rmt(T), Color,k)axis(0 10*T -200 200)set(gca,XTick,0:T:10*T) %给定x轴的刻度set(gca,XTickLabel,0,1,2,3,4,5,6,7,8,9,10) %给定x轴每个刻度的标示值ylabel (fontsize14rmv(cm/s), Color,k)legend(拟合曲线1, 拟合曲线2, 数值计算曲线)grid on%end这时我们得到非常吻合的曲线图像3-3,此时，b=0.00068。这说明，在磁阻尼作用情况下与空气阻力情况相类似，我们只考虑阻尼系数与速度的一次方的关系是不够的，只有当讨论阻尼系数与速度二次方的关系时，才能得到更精确的和b，得出更加吻合的图像。function Q=funxz(t,x)T