物理学阻尼振动中阻尼系数的研究 大学毕业设计.docx
《物理学阻尼振动中阻尼系数的研究 大学毕业设计.docx》由会员分享，可在线阅读，更多相关《物理学阻尼振动中阻尼系数的研究 大学毕业设计.docx（36页珍藏版）》请在冰点文库上搜索。
物理学阻尼振动中阻尼系数的研究大学毕业设计
阻尼振动中阻尼系数的研究
摘要：
在我们的力学教材和其他的有关力学知识的书中，对于阻尼振动中阻尼系数的研究已屡见不鲜，并且有好多学者和老师都对它进行了分析和研究，并取得了不错的理论结果。
但在我们课本中只是利用了传统的试验方法来测定阻尼振动中的阻尼系数，在本文中，我们实验测得的数据进行拟合来计算空气阻尼和磁阻尼的阻尼系数，得出了令人满意的结果，与传统实验手段比较起来更为准确。
关键词：
阻尼系数；阻尼振动；空气阻尼；磁阻尼；曲线拟合
Astudyondampingcoefficientofvibrationdamping
ABSTRACT:
Inourmechanicalmaterialsandotherbooksrelatingtothemechanicalknowledge,fordampingvibrationdampingconstantshavebeenfrequent,andtherearemanyscholarsandteachersofhiswereanalyzedandstudiedandmadegoodtheoreticalresults.Butinourtextbooks,justusethetraditionaltestmethodstomeasurethedampinginvibrationdamping,inthispaper,weusecurvefittingwiththeexperimentaldatacombinedwiththespecificproceduresinMatlabtocalculatetheairdampingandmagneticdamping,thedampingcoefficient,obtainedsatisfactoryresults,andexperimentalmeansalotmoreconvenient.
KEYWORDS:
Dampingcoefficient;dampingvibration;airdamping;magneticdamping;curvefitting
引言
机械振动是力学中的常见现象，实际生活中的振动并不是像理想条件下的简谐振动，而是带有阻尼的阻尼振动。
最近研究阻尼对振动的影响也成为人们比较关注的问题。
在实验上传统的测定方法是通过利用半衰期或利用阻尼振动曲线求出阻尼振子的阻尼系数
。
在本文中我们根据实验数据利用Matlab做出理论曲线并同时画出速度峰值的拟合曲线，再将二者进行比较估计的方法求得阻尼系数
。
在本文中，我们不仅利用这种方法讨论只在空气阻力下的阻尼系数并且画出拟合曲线，而且我们还讨论了在磁阻尼作用下的阻尼系数。
文中不仅讨论了阻尼系数在两种情况下与速度一次方的关系，还讨论了与速度二次方的关系。
测量过程简单易行且得到的实验数据可靠,具有一定的推广价值。
在下面我们分别对仅存在空气阻尼及存在磁阻尼的情况进行讨论。
1数据的曲线拟合
工程技术实践中,特别是在实验测试中,往往只能得到两个相关变量的一系列离散值,它们间的函数关系就只得用列表法或图示法表示。
但是，有时因为某种原因希望把这种关系转换成解析表达式。
另外,有些变量间的关系虽然可以用解析法表示，由于数学解析式过于繁杂不便使用，也希望用一个简单的解析表达式来近似地代替它，凡此种种，利用简单解析式
近似地代替列表法、图示法或复杂解析式表示的函数
一类问题，即寻找一个满足
的简单函数
的问题，都可称之为函数逼近。
数据拟合法就是函数逼近的重要方法之一。
科学实验或数据统计中，在得不到两组相关数据之间的精确解析表达式时，就希望找到关系
的近似解析式
，
越简单越好。
在这种情况下，我们提出了函数逼近方法-----数据拟合法，它不要求构造的近似函数
全部通过样本点，而是“很好的逼近”它们。
这种近似函数
反映了已知数据组间存在着某种关系的一般趋势，是“拟合”这些数据得出的函数曲线，不同于“插值”得出的曲线.寻找“很好的逼近”的函数有多种方法，我们常用的是最小二乘法。
在MatLab中，有专用的拟合指令，指令是polyfit，使用格式为：
p=polyfit（x,y,m），可以方便地给出求算结果，我们便利用实验得到的振动数据结果利用拟合方法，对数据点进行拟合，得到拟合曲线。
2模型
对于存在阻尼力的情况下，阻尼振动[1]方程为：
（2.1）
其中,
、
分别为系统的阻尼系数和固有频率.方程（2.1）的解为：
（2.2）
其中，
、
分别为弹簧振子的振幅、初始相位角，而
，
为阻尼振动周期
（2.3）
由方程（2.2）可得其速度与时间的关系式为：
（2.4）
这为我们的实验测量[2]提供了理论依据，根据公式所有公式都与
有关，因此只需要记录时间，周期即可。
3无磁阻尼时空气阻尼的阻尼系数的研究
实验所采用的气垫导轨型号是L-QG-Y-1500/5.8，电子计数器型号是MUJ-5C。
在气垫导轨上测量弹簧振子振动时的阻尼系数，进而求得时间常数
及品质因数
等物理量。
假定振动为线性振动，且气垫导轨存在的空气阻尼很小。
测量时所用滑块质量为179.70g，振动的振幅为30cm。
由实验测得弹簧振子振动连续10个周期的数据，见表3-1。
表3-1振动周期测量数据
n
1
2
3
4
5
6
7
8
9
10
T（s）
1.2690
1.2691
1.2691
1.2691
1.2690
1.2689
1.2690
1.2688
1.2687
1.2687
利用Matlab计算10个周期测量值的算术平均值、实验平均值的标准误差。
计算程序如下：
>>clear
hn=[1.2690,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;
forn=1:
N;
hp=hp+hn（n）/N;
end
hp
forn=1:
N;
vhn（n）=abs（hn（n）-hp）^2;
end
sigma2=0;
forn=1:
N;
sigma2=sigma2+vhn（n）/（N-1）;
end
sigman=sqrt（sigma2）;
s=sigman/sqrt（N）;
hp,s
%end
hp代表算术平均值,s为平均值的实验标准误差。
输出结果：
hp=1.2689
s=4.9889e-005
即：
.
在实验时取初始位移
，初始速度为
，并通过放置在平衡位置的电子计数器记录下连续10个周期中，一个周期从正方向
和反方向
通过两次平衡位置时光电门的速度。
表3-2滑块通过平衡位置时的速度
0.25
1.25
2.25
3.25
4.25
5.25
6.25
7.25
8.25
9.25
145.98
142.45
139.86
137.36
134.95
132.28
129.87
127.71
125.63
123.30
0.75
1.75
2.75
3.75
4.75
5.75
6.75
7.75
8.75
9.75
-143.88
-141.04
-138.50
-136.05
-133.69
-131.06
-128.86
-126.74
-124.53
-122.25
程序设计：
clear
T=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）;%生成拟合多项式（向量形式）,阶次为3
pvt1=poly2str（pv,'t'）;%生成拟合多项式（字符形式）
pvt1
输出结果：
pvt1=
0.032136t^2-2.3349t+146.4178
clear
T=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）;%生成拟合多项式（向量形式）,阶次为2
pvt2=poly2str（pv,'t'）;%生成拟合多项式（字符形式）
pvt2
输出结果：
pvt2=
-0.021055t^2+2.1588t-145.8269
，
分别是沿正方向运动时，通过平衡位置的光电门的时间和速度。
，
是沿负方向运动时，通过平衡位置的光电门的时间和速度。
3.1阻尼系数
由阻尼振动方程的解（2.4）式编制Matlab程序计算理论上时间与速度曲线,同时将两拟合曲线绘制在同一幅图中。
程序设计如下：
>>clear
A=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=-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（'\fontsize{14}\rmt（T）','Color','k'）
axis（[010*T-200200]）
set（gca,'XTick',0:
T:
10*T）%给定x轴的刻度
set（gca,'XTickLabel',{'0','1','2','3','4','5','6','7','8','9','10'}）%给定x轴每个刻度的标示值
ylabel（'\fontsize{14}\rmv（cm/s）','Color','k'）
legend（'拟合曲线1','拟合曲线2','理论曲线'）
gridon
%end
该程序绘制出弹簧振子时间与速度曲线,见图2-1与2-2。
其中实线、点线是由实验上得到向正方向、负方向通过平衡位置光电门的时间与速度数据拟合的曲线,实线是由理论方程（2.3）所得到的曲线[3]。
即为所要求的阻尼系数，在用Matlab画图时
的值是在不断地调试的，当将解中的
调整为
时,理论上曲线通过平衡位置时的速度与两拟合曲线基本吻合,即两按拟合曲线正是速度曲线的包络线。
所以得到此时阻尼系数为：
（4）
可以看出通过对
的多次取值估计，重复多次计算画图，我们得出了较为满意阻尼系数，使得理论线与拟合曲线基本吻合，如图2-1所示，结果证明上述理论与实验结合的方法是可行的。
与此同时，我们还可以得到时间常数和品质因数：
时间常数为
（5）
品质因数为
（6）
但是我们从图2-2看出拟合曲线并没有很好的与理论实线相吻合，放大后观察到还有不相交的部分，也就是不完全包络，原因在于我们绘制的拟合曲线是一次直线，所以结果并不是令人很满意。
；
无论对于
怎么调整，对于放大的图像来看，总不可能全部包络速度的峰值，只能是尽可能的平衡，这说明仅仅调整
是不够的，还要考虑阻尼与速度的二次方的关系。
3.2考虑阻尼与速度二次方有关
上面所讨论的是空气阻尼与速度一次方的拟合关系，在某些点上并不能很好的拟合，下面利用MatLab来模拟与速度二次方[5]的拟合曲线。
在M文档中建立如下函数：
functionQ=funxz（t,x）
T=1.2689;
omegaf=2*pi/T;
d=0.0145;
p=0.00043;
Q=[x
（2）;-omegaf^2*x
（1）-2*d*x
（2）-2*p*x
（2）.^2];
程序如下：
>>clear
[t,x]=ode45（@funxz,[0:
0.01:
14],[-300]）;
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','LineWidth',2）
xlabel（'\fontsize{14}\rmt（T）','Color','k'）
axis（[010*T-200200]）
set（gca,'XTick',0:
T:
10*T）%给定x轴的刻度
set（gca,'XTickLabel',{'0','1','2','3','4','5','6','7','8','9','10'}）%给定x轴每个刻度的标示值
ylabel（'\fontsize{14}\rmv（cm/s）','Color','k'）
legend（'拟合曲线1','拟合曲线2','数值计算曲线'）
gridon
%end
在M文档中所建立的函数恰好是阻尼系数与速度二次方的函数关系式，我们又对
和p进行了多次调试，结果如图2-3与图2-4所示，通过对阻尼与速度二次方关系的讨论，所得结果分别与图2-1和图2-2比较，拟合曲线与速度峰值更好的吻合了，并且得到了更为精确的
，并且得出了速度二次方因子p=0.00043。
；
。
4研究磁阻尼情况下的阻尼系数
在滑块两侧对称地粘贴上永久小磁铁块，此时由于磁阻尼[4]力的作用，使振动的振幅衰减的更快。
在此时测量的阻尼系数，仍假定振动为线性振动。
测量时所用滑块质量为191.06g，振动的振幅为30.00cm。
由实验测得弹簧振子连续振动10个周期的数据，见表4-1。
表4-1振动周期测量数据
n
1
2
3
4
5
6
7
8
9
10
T（s）
1.3000
1.3000
1.2996
1.2998
1.2996
1.2999
1.2999
1.2998
1.2997
1.2998
利用MatLab计算10个周期测量值的算术平均值、实验平均值的标准误差。
计算程序：
>>clear
hn=[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;
forn=1:
N;
hp=hp+hn（n）/N;
end
forn=1:
N;
vhn（n）=abs（hn（n）-hp）;
end
forn=1:
N;
vhn2（n）=（hn（n）-hp）^2;
end
sigma2=0;
forn=1:
N;
sigma2=sigma2+vhn2（n）/（N-1）;
end
sigman=sqrt（sigma2）;
s=sigman/sqrt（N）;
hp,s%hp代表算术平均值,s为平均值的实验标准误差.
%end
hp=1.2998
s=4.5826e-005
即：
。
在实验时取初始位移
，初始速度为
。
我们启记录下连续10个周期中，一个周期从正方向和反方向通过两次平衡位置时光电门的速度。
表4-2振动周期测量数据
0.25
1.25
2.25
3.25
4.25
5.25
6.25
7.25
8.25
9.25
140.65
135.32
131.06
126.74
122.85
119.05
114.94
111.23
107.76
104.49
0.75
1.75
2.75
3.75
4.75
5.75
6.75
7.75
8.75
9.75
-137.74
-133.33
-128.86
-124.84
-120.92
-117.23
-113.25
-109.41
-106.04
-102.77
，
是向正方向运动时，通过平衡位置的光电门的时间和速度。
根据上述数据编程拟合出
，
的关系式。
程序设计：
>>clear
T=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,104.49];%输入速度
pv=polyfit（te,ve,2）;%生成拟合多项式（向量形式）,阶次为3
pvt1=poly2str（pv,'t'）;%生成拟合多项式（字符形式）
pvt1
输出结果：
；
由上述结果可得
与
之间的关系为：
、
是向负方向运动时，通过平衡位置的光电门的时间和速度。
程序设计：
>>clear
T=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.41,-106.04,-102.77];%输入速度
pv=polyfit（te,ve,2）;%生成拟合多项式（向量形式）,阶次为3
pvt2=poly2str（pv,'t'）;%生成拟合多项式（字符形式）
pvt2
输出结果：
。
由上述结果可得
与
之间的关系为：
4.1阻尼系数
由阻尼振动方程（3）式的解编制Matlab程序计算理论上给出的时间与速度曲线，同时将两拟合曲线绘制在同一幅图中。
程序中的vn1即
，vn2即
。
程序设计：
>>clear
A=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+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（'\fontsize{14}\rmt（T）','Color','k'）
axis（[010*T-200200]）
set（gca,'XTick',0:
T:
10*T）%给定x轴的刻度
set（gca,'XTickLabel',{'0','1','2','3','4','5','6','7','8','9','10'}）%给定x轴每个刻度的标示值
ylabel（'\fontsize{14}\rmv（cm/s）','Color','k'）
legend（'拟合曲线1','拟合曲线2','理论曲线'）
gridon
%end
该程序绘制出弹簧振子时间与速度曲线，见图3-1。
其中点划线、点线是由实验上得到向正方向、负方向通过平衡位置光电门的时间与速度数据拟合的曲线，实线是理论方程所得到的曲线。
当将解中的
调整为
时，理论上曲线通过平衡位置时的速度与两拟合曲线基本吻合。
由此得到阻尼系数为
（4）
时间常数
为
（5）
品质因数
为
（6）
；
如图3-2.、图3-1所示，我们得到了基本吻合的拟合曲线，但放大后观察，不论
如何调整，还是没有得到满意的效果，那么就要考虑与速度二次方的关系。
4.2考虑速度二次方因子
上面所讨论的是磁阻尼与速度一次方的拟合关系，在某些点上并不能很好的拟合，下面利用MatLab来模拟与速度二次方[5]的拟合曲线。
在M文档中建立如下函数：
functionQ=funxz（t,x）
T=1.2998;
omegaf=2*pi/T;
d=0.025;
b=0.00068;
Q=[x
（2）;-omegaf^2*x
（1）-2*d*x
（2）-2*b*x
（2）.^2];
程序如下：
>>clear
[t,x]=ode45（@funxz,[0:
0.01:
14],[-300]）;%[t,x]=ode45（@funerx,[0:
0.01:
2*pi],[40]）;
T=1.2998;
t1=0:
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（'\fontsize{14}\rmt（T）','Color','k'）
axis（[010*T-200200]）
set（gca,'XTick',0:
T:
10*T）%给定x轴的刻度
set（gca,'XTickLabel',{'0','1','2','3','4','5','6','7','8','9','10'}）%给定x轴每个刻度的标示值
ylabel（'\fontsize{14}\rmv（cm/s）','Color','k'）
legend（'拟合曲线1','拟合曲线2','数值计算曲线'）
gridon
%end
这时我们得到非常吻合的曲线图像3-3,此时
，b=0.00068。
这说明，在磁阻尼作用情况下与空气阻力情况相类似，我们只考虑阻尼系数与速度的一次方的关系是不够的，只有当讨论阻尼系数与速度二次方的关系时，才能得到更精确的
和b，得出更加吻合的图像。
functionQ=funxz（t,x）
T
升级会员 