单自由度系统自由振动matlab程序

山东大学Ｍatｌａb 课程作业学院：机械工程学院专业：姓名:学号:基于Simｕｌink仿真得振动学问题解决实例1.单自由度无阻尼自由振动仿真表达式：仿真框图：参数设置：k=1０0N/m m＝4kｇ初始状态：初速度为0 初始位移为５仿真结果：2.简谐波形得里沙茹图形分析仿真框图：参数设置：K=100m＝４→raｄ／ｓSin wａve参数设置：Ａｍplｉtｕde1 ；Fｒeqｕencｙ 5 １０15初始状态：①→φ=②→φ=③＝１，=5→φ=45；④＝1，＝−5→φ=135;⑤＝0,=−1→φ=１80ＸY Graph参数x-mｉn －2；x-max 2；y－min—2; y-ｍａx 2Frｅquency 5时仿真结果:Ｆrｅquenｃy 1０时仿真结果：Ｆreｑｕency 15时仿真结果：3.单自由度有阻尼自由振动表达式:仿真框图：参数设置：ﻫ令k=１00,m＝１0，ｃ=10 初始状态：ﻫ初始速度为0，位移为1仿真结果：４、衰减振荡得阻尼比得估计参数：ｋ=1００,m=10，c=２初始条件:x0=1，v0＝0仿真图框：初始振幅为1，约7个周期时衰减为０、25,对数减幅：δ＝（ｌn４）/７≈0、０99阻尼比§≈δ/2≈0、0３2理论值§＝0、5c(km)−0、5≈0、０325、单自由度有阻尼+正弦激励表达式：令激励则方程变形为参数设置：令k=４，ｍ=1，c=０、2初始状态：ﻫ初始速度为0,位移为0、０5 仿真框图：仿真结果：6、利用速度共振得里沙茹图进行固有频率与阻尼系数分析仿真框图:改变激励频率：=1、２；1、6;1、8；１、9;１、95；2；2、05;2、1;2、2等7、两自由度无阻尼系统自由振动表达式:参数设置:m1=１，m2=2 k１=1,ｋ2=1，k３=2初始状态:①速度0，ｍ1、m2位移均为１②速度0，m１位移1，m2位移−0、5③速度0,m1位移1，m2位移0 仿真结果：①②③。

采用EL-CENTRO地震波计算单自由体系的振动位移反应谱前言：本课程论文采用软件是MATLAB R2010a，所编程序应用的方法是线性加速度法。

程序采用的是两层循环上实现位移谱的求取。

内循环实现在地震波下由上一步所得到的位移、速度和加速度得到下一时间间隔后的位移速度和加速度；外循环实现在不同结构自振周期（频率）下遭遇地震波的响应。

本文在求出相对位移反应谱的同时，也求出了相对速度反应谱、绝对加速度反应谱。

一．线性加速度法的简述：线性加速度法是直接数值积分法求解地震反应的方法之一，本文所采用的线性加速度法参考大崎顺彦的《地震动的谱分析入门》第二版。

具体计算公式详见大崎顺彦的《地震动的谱分析入门》第二版P116-P118。

二．所编程序及编译：clear% ***********读入地震记录***********fid = fopen('ei.txt');[Accelerate,count] = fscanf(fid,'%g'); %count 读入的记录的量time=0:0.02:(count-1)*0.02;% ***********线性加速度法计算各反应***********%初始化各储存向量Displace=zeros(1,count); %相对位移Velocity=zeros(1,count); %相对速度AbsAcce=zeros(1,count);%绝对加速度Damp=0.05; %结构阻尼比取为0.05Tc=0.0:0.05:10; %结构自振周期Dt=0.02; %地震记录的步长%记录计算得到的反应，MDis为最大相对位移，MVel为最大相对速度%MAcc某阻尼时最大绝对加速度，用于画图MDis=zeros(1,length(Tc));MVel=zeros(1,length(Tc));MAcc=zeros(1,length(Tc));t=1; %在下一个循环中控制不同的结构自振周期for T=0.0:0.05:10Frcy=2*pi/T ; %结构自振频率DamFrcy=Frcy*sqrt(1-Damp*Damp);%计算公式化简e_t=exp(-Damp*Frcy*Dt);s=sin(DamFrcy*Dt);c=cos(DamFrcy*Dt);A=zeros(2,2);A(1,1)=e_t*(s*Damp/sqrt(1-Damp*Damp)+c);A(1,2)=e_t*s/DamFrcy;A(2,1)=-Frcy*e_t*s/sqrt(1-Damp*Damp);A(2,2)=e_t*(-s*Damp/sqrt(1-Damp*Damp)+c);d_f=(2*Damp^2-1)/(Frcy^2*Dt); %计算公式化简d_3t=Damp/(Frcy^3*Dt);B=zeros(2,2);B(1,1)=e_t*((d_f+Damp/Frcy)*s/DamFrcy+(2*d_3t+1/Frcy^2)*c)-2*d_3t;B(1,2)=-e_t*(d_f*s/DamFrcy+2*d_3t*c)-1/Frcy^2+2*d_3t;B(2,1)=-e_t*(((Damp/(Frcy*Dt)+1)*s/DamFrcy)+(1/(Frcy^2*Dt))*c)+1/(Frcy^2*Dt);B(2,2)=e_t*((Damp/(Frcy*Dt)*s/DamFrcy)+(1/(Frcy^2*Dt))*c)-1/(Frcy^2*Dt);for i=1:(count-1) %根据地震记录,计算不同的反应Displace(i+1)=A(1,1)*Displace(i)+A(1,2)*Velocity(i)+B(1,1)*Accelerate(i)+B(1,2)*Accelerate(i +1);Velocity(i+1)=A(2,1)*Displace(i)+A(2,2)*Velocity(i)+B(2,1)*Accelerate(i)+B(2,2)*Accelerate(i+ 1);AbsAcce(i+1)=-2*Damp*Frcy*Velocity(i+1)-Frcy^2*Displace(i+1);endMDis(1,t)=max(abs(Displace));MVel(1,t)=max(abs(Velocity));if T==0.0MAcc(1,t)=max(abs(Accelerate)); %当结构的自振周期为0时，其绝对加速度应等于地面加速度elseMAcc(1,t)=max(abs(AbsAcce));endDisplace=zeros(1,count);%初始化各储存向量避免下次计算时引用到前面结果Velocity=zeros(1,count);AbsAcce=zeros(1,count);t=t+1; %t=length(Tc)，即所求结构自振周期有多少个，对应就运行多少次。

振动力学是研究物体在作往复振动或周期性运动时的力学规律和特性的一门学科。

它在工程、物理、地震学等领域中有着广泛的应用。

MATLAB是一种强大的数值计算和科学绘图软件，可以用于振动力学的建模、仿真和可视化。

在振动力学基础方面，需要掌握以下内容：
1. 单自由度系统：这是振动力学的基础，主要研究质点的简谐振动和阻尼振动等。

需要了解自由度、刚度、阻尼和质量等概念，并能够利用牛顿第二定律、欧拉-拉格朗日原理等方法分析运动方程和相应的振动特性。

2. 多自由度系统：多自由度系统是复杂振动问题的常见形式，需要掌握刚体系统、弹性系统和连续系统等的振动特性。

这里需要了解模态分析、正交性原理和频率响应等概念，并学会通过欧拉-拉格朗日方程和质量矩阵、刚度矩阵等进行系统参数的求解和模拟。

在MATLAB应用方面，需要掌握以下内容：
1. MATLAB基础语法和常用命令，如数据类型、矩阵运算、函数定义和图形绘制等。

2. 振动力学的MATLAB模型建立和仿真分析。

需要学会利用MATLAB解决振动力学问题的程序设计和编写，如求解ODE方程组、进行模态分析和频率响应分析等。

3. MATLAB可视化工具的使用，如画图工具箱、动画工具箱、GUI界面设计与应用等，以便更加直观地展现振动力学问题的结果和结论。

振动力学基础与MATLAB应用是一门需要深入掌握的学科。

通过深入学习这门学科，可以更好地理解和应用振动力学的理论和方法，同时也可以更好地掌握MATLAB在振动力学中的应用。

1.function f = hanning_imp(t, Tc, A)2. f = zeros(size(t));3. f(t < Tc) = A / 2 * (1 - cos(2*pi * t(t < Tc) / Tc));4.end下面是计算单自由度系统响应的Matlab程序，计算传递函数，画实频、虚频、幅频、相位、导纳图[plain]view plain copy1.m = 100;2.k = 1000;3. c = 100;4.5.num = 1;6.den = [m c k];7.sys = tf(num, den);8.dt = 0.00001;9.fs = 1/dt; %采样频率(Hz) 100Hz 实际并不需要这么高的采样频率，但是如果采样时间太小，hanning脉冲不完整10.% 为了得到准确的响应dt一定要小，否则做出的相位可能不对11.12.t = 0:dt:200;13.Tc = 0.001;14. A = 10;15.u = hanning_imp(t, Tc, A);16.y = lsim(sys, u, t);17.% y = impulse(sys, t);18.y = y';19.20.N = length(u);21.fy = fft(y);22.fu = fft(u);23.ft = fy ./ fu;24. f = (0:N-1) * fs ./ N;25.ft_r = real(ft);26.ft_i = imag(ft);27.28.part = (f < 30);29.30.figure('name', '实频');31.plot(f(part), ft_r(part));32.figure('name', '虚频');33.plot(f(part), ft_i(part));34.figure('name', '幅值');35.plot(f(part), abs(ft(part)));36.figure('name', '相位');37.plot(f(part), phase(ft(part)) * 180 / pi);38.figure('name', '导纳圆');39.plot(ft_r(part), ft_i(part), '.');40.axis equal;41.hold on;42.xk = [-0.000116543. 0.000285744. 8.297e-545. 0.000486946. 0.000682];47.xk1 = [-0.000801748. -0.000308249. 0.000285750. 0.00119451. 0.001451];52.yk1 = [-0.00281753. -0.00309554. -0.00318555. -0.00286956. -0.002633];57.58.yk = [-0.00314959. -0.00318560. -0.00317961. -0.00316562. -0.003121];63.k = 5;64. A = [65. sum(xk.^2) sum(xk .* yk) sum(xk)66. sum(xk.*yk) sum(yk.^2) sum(yk)67. sum(xk) sum(yk) k];68. B = -[69. sum(xk.^3 + xk.* yk.^2)70. sum(xk.^2.*yk + yk.^3)71. sum(xk.^2+yk.^2)];72.rlt = A\B;73.x0 = rlt(1)*(-0.5);74.y0 = rlt(2)*(-0.5);75.r = sqrt(rlt(1)^2/4+rlt(2)^2/4-rlt(3));76.fai = 0:0.01:2*pi;77.x = x0 + r * cos(fai);78.y = y0 + r * sin(fai);79.plot(x,y, 'r');实验模态分析------非数学公式的简单概述之二分类：模态空间译文| 标签：锤击法激振器窗函数曲线拟合2011-07-02 23:00阅读(3040)评论(0)为何只需获得频响函数矩阵的一行或一列？理解从可能得到频响函数矩阵的不同元素中得到模态振型对我们来说是非常重要的。

funct ion responsenew%基于戏线性滞冋模型的单自由度体系的地震能虽分析程序%made by Yanan Li Faculty of Infrastructure Enginoering弔Dalian University of Technology%质量57041kg,阻尼36612 N s/叫初始刚度2350000N/m,刚度折减系数0. 2,屈服位移0・01m, 采用ELCENTRO波%参数替换直接在下面修改，然后运行clcformat long;m=57041;% 质量ug= importdataC* ELCENTRO. txt') ;%地震波txt 文件ug=ug/100;P=-m*ug;num=size(P, 1);c 二36612; % 阻尼k 1=2350000; %初始刚度k2二kl*0・2;a=zeros(l t num) ;v=zeros(l, num) ;x=zeros(l, num);%加速度速度位移Ei=zeros(l, num);Ek二zeros(1, num);Ed二zeros(1, num); Eh二zeros (1, num)输入能动能阻尼耗能滞回耗能EI=zeros(l t num); EK二zeros (1, num) ;ED=zeros(l, num) ;EH=zeros(l, num) ;%累积的各种能量t imc=zeros(l, num);a(l)=P (l)/m;hfl=zeros (1, num);h=0. 02;%地震波采样间隔Xy=0.01;%屈服位移pxmax=0;nxmax-0;pd=0;%双线型滞冋模型折线的标识，0农示弹性，1表示正向弹犁性，2表示反向弹性，-1 表示反向弹塑性，-2表示正向弹性for ii=l:numif pd==0 %弹性阶段k二kl;if x(ii)>Xypd=l;b=[(a(ii)-a(ii-l))/6/h a(ii~l)/2 v(ii-l) x(iiT)-Xy];% 拐点处理d=roots(b);for j=l:length(d)if isreal (d(j))==ldth=d(j);endendhp二h-dth;vp=v(ii-l)+a(ii-l)*dth+((a(ii)-a(ii-l))*dth*2/h/2);ap=a(i i-l) + (a (i i)-a(ii-l)) *dth/h;pp=m*ap+c*vp+kl*Xy;kd=k2+3*c/hp+6*m/hp/hp;dtpd=P(i i +1) -pp+m* (6*vp/hp+3*ap) +c* (3*vp^-hp*ap/2); dtx=dtpd/kd;dtv=3*dtx/hp-3*vp-hp*ap/2;x(ii)=Xy+dtx;v(ii)=vp^dtv;dtf=k2*dtx;hfl(ii)=kl*Xy+dtf;a(ii) = (P(ii)-hfl (i i+1)-c*v(i i))/m;elseif x(ii)<-Xypd=-l;b=[(a(ii)-a(ii-l))/6/h a(ii-l)/2 v(iiT) x(iiT)+Xy];% 拐点处理d=roots(b);for j=l:length(d)if isreal (d(j))==ldth=d(j);endendhp=h-dth;vp=v(i i 1)・a(i i 1)*dth+ ((a(ii)-a(ii-l))*dth"2/h/2); ap a(ii-l) + (a(ii)-a(ii-l))*dth/h;pp=m*ap+c*vp_k1*Xy; kd=k2+3*c/hp+6*m/hp/hp;dtpd=P (i i + 1)-pp+m*(6*vp/hp+3*ap)(3*vp^hp*ap/2);dtx=dtpd/kd;dtv-3*dtx/hp-3*vp-hp*ap/2;x(ii)=-Xy+dtx： v(i i)=vp^dtv;dtf=k2*dtx; hfl(ii)=-kl*Xy+dtf;a(i i) = (P(i i)-hf 1 (i i)-c*v(ii))/m;endendif pd==l %正向弹塑性k=k2;if v(ii)<0pd=2;b=[(a(ii)-a(ii-l))/2/h a(ii-l) v(ii-l)];% 拐点处理d二rools(b);for j=l:length(d)if isreal (d(j))==ldth=d (j);endendhp=h-dth;xp=x(iiT)+v(i i-1) *dth+a (i i-1) *dth*2/2+ (a (i i) -a (i i-l))*dth*3/h/6; ap=a (i i-1)+ (a (i i) -a (i i-1)) *dth/h;pp-m*ap+c*O+k1*Xy+k2* (xp-Xy);kd=kl+3*c/hp+6*m/hp/hp:dtpd=P(ii+l)-pp+m*(6*0/hp+3*ap)+c*(3*0+hp*ap/2);dtx=dtpd/kd;dtv=3*dtx/hp-3*0-hp*ap/2;x(i i)=xp+dtx;v(ii)=O+dtv;dtf=kl*dtx;hfl(ii)=kl*Xy+k2*(xp-Xy)+dtf;a(ii) = (P(ii)-hfl (ii)-c*v(ii))/m;pxmax=xp;endif pd==2 %反向弹性k=kl;if x (i i)>pxmaxpd=l; ■b=[(a(ii)-a(ii-l))/6/h a(ii-l)/2 v(ii-l) x(ii-l)-pxmax]:% 拐点处理d=roots(b);for j=l:length(d)if isreal (d(j))==I dth=d(j);endendhp=h-dth; vp=v(ii-l)+a(ii-l)*dth+((a(ii)-a(ii-l))*dth"2/h/2); ap=a(ii-l)+(a(ii)-a(ii-l))*dth/h;pp-m*ap+c*vp+k1*Xy+k2*(pxmax-Xy); kd=k2+3*c/hp+6*m/hp/hp;dtpd=P(ii+1)-pp+m*(6*vp/hp+3*ap)+c*(3*vp+hp*ap/2); dtx=dtpd/kd;d t v=3*d1x/hp-3*vp-hp*ap/2;x(i i)二pxmax+dtx;v(ii)=vp+dtv;dtf=k2*dtx;hfl (i i)=kl*Xy+k2*(pxmax-Xy) +dtf;a(ii) = (P(ii)-hfl (ii)-c*v(ii))/m;el seif x (i i)<(pxmax~2*Xy)pd=-l;―.・pu og 、(( =)A *3—(二二Jq —(二)d )H(二)P,±P +AX盖〒Ax*上—殳艾舌 X F (二二jqEp*znp 2p+dA:=(=)A-xlp+AX*zlxpulxdu(二)x,z 、dB*dqld 家gldv、xlp^HA4p&<p d "P H X 4p-(z 、dp*dfdA*g)芯+(du*g+dq、dA*9)*u+dd —(7二)dupdpp -dq'dq'lu 关 9+dqb苗+0上HP上-(xx*^IAX*>l —z>l*xeulxd) +d A *o +d e *u "dd2>-p*(〒二)E —(=)B)+U丄一)fde二z 'q 'Fe p 決((TT OF(U)P ))+±P*II )P +UI )AH D A召p —qudqp u a)pUO」(「)p*pTH((「)P 二p冒岛(P)£8U 上二丄joj-(q)s_oof p三P+O A=)Awp+dxH(二)x-z 、dydq —0^—dq、w'pdup,(z 、d¥dq+o為)*¥(dp*e+dq/0*9)*lu+dd —(I +二)dupdlp‘ d<dq 、ul*9+d<。

一、引言在工程领域中，单自由度振动系统是一种常见的动力学模型，其在建筑结构、机械设备等领域都有重要的应用。

而计算单自由度振动系统的反应是工程设计和分析中的重要任务之一，对于系统的稳定性、安全性等具有重要意义。

为了快速而准确地计算单自由度振动系统的反应，工程师和研究人员常常使用MATLAB编写程序来实现计算。

二、MATLAB在单自由度振动反应计算中的优势1. 灵活性：MATLAB是一种功能强大的编程语言和工具，可以实现复杂的算法和数学模型，能够满足工程设计和分析中的各种需求。

2. 可视化：MATLAB具有丰富的绘图和可视化功能，可以直观地展示单自由度振动系统的反应结果，使工程师和研究人员更好地理解系统的运动特性。

3. 高效性：MATLAB提供了丰富的计算和求解工具，可以快速而准确地计算单自由度振动系统的反应，节约了工程师和研究人员的时间和精力。

三、MATLAB编写单自由度振动反应计算程序的基本步骤1. 确定系统参数：首先需要确定单自由度振动系统的质量、刚度、阻尼系数等参数，这些参数将影响系统的振动响应。

2. 构建系统模型：根据系统的参数和运动方程，可以利用MATLAB编写对应的单自由度振动系统模型。

3. 求解运动方程：利用MATLAB提供的求解工具，可以求解单自由度振动系统的运动方程，得到系统的振动响应。

4. 可视化结果：最后可以利用MATLAB的绘图和可视化工具，将系统的振动响应以图表的形式展现出来，便于工程师和研究人员对系统的运动特性进行分析和评估。

四、MATLAB编写单自由度振动反应计算程序的示例代码以下是一个简单的MATLAB示例代码，用于计算单自由度振动系统的阻尼比为0.2时的阻尼比对应的阻尼比比，供读者参考。

```matlab定义系统参数m = 1; 质量k = 10; 刚度zeta = 0.2; 阻尼比omega_n = sqrt(k / m); 自然频率计算阻尼比对应的阻尼比比omega_d = omega_n * sqrt(1 - zeta^2);disp(['阻尼比对应的阻尼比比为：', num2str(omega_d /omega_n)]);```通过上述示例代码，可以看出MATLAB的编写方式简单明了，利用MATLAB可以快速计算单自由度振动系统的各种响应参数。

1.function f = hanning_imp(t, Tc, A)2. f = zeros(size(t));3. f(t < Tc) = A / 2 * (1 - cos(2*pi * t(t < Tc) / Tc));4.end下面是计算单自由度系统响应的Matlab程序，计算传递函数，画实频、虚频、幅频、相位、导纳图[plain]view plain copy1.m = 100;2.k = 1000;3. c = 100;4.5.num = 1;6.den = [m c k];7.sys = tf(num, den);8.dt = 0.00001;9.fs = 1/dt; %采样频率(Hz) 100Hz 实际并不需要这么高的采样频率，但是如果采样时间太小，hanning脉冲不完整10.% 为了得到准确的响应dt一定要小，否则做出的相位可能不对11.12.t = 0:dt:200;13.Tc = 0.001;14.A = 10;15.u = hanning_imp(t, Tc, A);16.y = lsim(sys, u, t);17.% y = impulse(sys, t);18.y = y';19.20.N = length(u);21.fy = fft(y);22.fu = fft(u);23.ft = fy ./ fu;24.f = (0:N-1) * fs ./ N;25.ft_r = real(ft);26.ft_i = imag(ft);27.28.part = (f < 30);29.30.figure('name', '实频');31.plot(f(part), ft_r(part));32.figure('name', '虚频');33.plot(f(part), ft_i(part));34.figure('name', '幅值');35.plot(f(part), abs(ft(part)));36.figure('name', '相位');37.plot(f(part), phase(ft(part)) * 180 / pi);38.figure('name', '导纳圆');39.plot(ft_r(part), ft_i(part), '.');40.axis equal;41.hold on;42.xk = [-0.000116543. 0.000285744. 8.297e-545. 0.000486946. 0.000682];47.xk1 = [-0.000801748. -0.000308249. 0.000285750. 0.00119451. 0.001451];52.yk1 = [-0.00281753. -0.00309554. -0.00318555. -0.00286956. -0.002633];57.58.yk = [-0.00314959. -0.00318560. -0.00317961. -0.00316562. -0.003121];63.k = 5;64.A = [65. sum(xk.^2) sum(xk .* yk) sum(xk)66. sum(xk.*yk) sum(yk.^2) sum(yk)67. sum(xk) sum(yk) k];68.B = -[69. sum(xk.^3 + xk.* yk.^2)70. sum(xk.^2.*yk + yk.^3)71. sum(xk.^2+yk.^2)];72.rlt = A\B;73.x0 = rlt(1)*(-0.5);74.y0 = rlt(2)*(-0.5);75.r = sqrt(rlt(1)^2/4+rlt(2)^2/4-rlt(3));76.fai = 0:0.01:2*pi;77.x = x0 + r * cos(fai);78.y = y0 + r * sin(fai);79.plot(x,y, 'r');实验模态分析------非数学公式的简单概述之二分类：模态空间译文| 标签：锤击法激振器窗函数曲线拟合2011-07-02 23:00阅读(3040)评论(0)为何只需获得频响函数矩阵的一行或一列？理解从可能得到频响函数矩阵的不同元素中得到模态振型对我们来说是非常重要的。

·kxcxckxmm xO 0cos F tw 实验名称 单自由度系统数值模拟一、实验目的、要求一、实验目的、要求1．熟悉单自由度系统强迫振动特性和求解方法；．熟悉单自由度系统强迫振动特性和求解方法； 2．掌握强迫振动系统的计算机模拟仿真方法。

．掌握强迫振动系统的计算机模拟仿真方法。

二、实验设备及仪器1. 计算机计算机2. Matlab 软件软件3. c 语言语言 三、实验步骤1．利用如右图所示的受力分析，得出单自由度系统强迫振动的运动方程。

自由度系统强迫振动的运动方程。

物体沿水平方向振动，取物体无扰力下的静平衡位置为坐标原点，水平向右为x 轴正向，建立如图所示的坐标系。

受力情况如图，其激励力为：0cos F F t w =，其中，，其中，0F 称为激励力的力幅，为常值。

称为激励力的力幅，为常值。

w 为激励频率，为常值。

为激励频率，为常值。

根据牛顿第二定律，得到单自由度系统强迫振动的运动方程：强迫振动的运动方程：0cos m x F t kx cxw =-- 2．对方程进行求解。

令n km w =，00F X k =，22c n c km m w ==，22c n c c c c m kmz w ===则原方程可以变形为：则原方程可以变形为：2202cos n n n x x x X t zw w w w ++=这是一个非齐次二阶常系数微分方程，根据微分方程理论，它的解由两部分组成这是一个非齐次二阶常系数微分方程，根据微分方程理论，它的解由两部分组成 12x x x =+其中，1x 代表齐次微分方程220n n x x x zw w ++=的解，简称齐次解，当1z <时，由前面的单自由度阻尼自由振动可得：前面的单自由度阻尼自由振动可得：()112cos sin cos()n n ttd d d x eBt B t Aet zwzww w w j --=+=-其中：21d n w zw =-×，称为衰减振动的固有频率。

Matlab技术在机械振动分析中的应用案例引言：机械振动作为机械工程领域中非常重要的研究方向，对于机械设备性能的评估和故障诊断具有关键作用。

随着计算机技术的飞速发展，Matlab作为一种功能强大的数学计算软件，被广泛应用于机械振动领域。

本文将通过介绍一些典型的应用案例，展示Matlab在机械振动分析中的优越性和实用性。

一、弹簧振动分析弹簧振动是机械系统中常见的一种振动形式。

通过Matlab可以方便地建立弹簧振动的数学模型，进行分析和仿真。

以弹簧单自由度系统为例，我们可以通过编写Matlab程序来求解该系统的振动特性，比如自然频率、阻尼比等参数。

此外，Matlab还提供了丰富的绘图功能，可以用来绘制系统的振动曲线和频谱图，进一步分析和评估系统的性能。

二、子午线摆振动分析子午线摆是一种简单而重要的振动系统，在物理实验教学中被广泛应用。

利用Matlab可以实现子午线摆的运动仿真和数据分析。

通过建立子午线摆的运动微分方程，我们可以利用Matlab的数值求解功能来模拟摆的运动过程，并绘制出摆角随时间的变化曲线。

此外，Matlab还可以计算出摆的周期和频率，提供了便捷的数据处理方法，方便进行实验数据的比对和验证。

三、转子系统振动分析转子系统的振动分析是机械工程中一项关键任务。

Matlab提供了大量的信号处理和频谱分析工具，可以用来对转子系统的动态性能进行评估和诊断。

首先，我们可以通过Matlab对转子系统的模态进行分析，求解出转子的模态频率和振型。

接着，利用Matlab的FFT函数进行频谱分析，可以得到转子系统的频谱图，并进一步分析出存在的谐波成分。

通过与参考频谱进行比较，我们可以判断转子系统是否存在异常振动，进而评估其工作状态。

四、车辆悬架系统振动分析车辆悬架系统的振动特性直接影响着驾驶员的驾驶感受和乘坐舒适度。

Matlab在车辆悬架系统的振动分析中发挥着重要作用。

通过建立车辆悬架系统的动力学模型，并利用Matlab进行模拟和仿真，我们可以得到车辆在不同路况下的振动响应。

funct ion responsenew%基于戏线性滞冋模型的单自由度体系的地震能虽分析程序%made by Yanan Li Faculty of Infrastructure Enginoering弔Dalian University of Technology%质量57041kg,阻尼36612 N s/叫初始刚度2350000N/m,刚度折减系数0. 2,屈服位移0・01m, 采用ELCENTRO波%参数替换直接在下面修改，然后运行clcformat long;m=57041;% 质量ug= importdataC* ELCENTRO. txt') ;%地震波txt 文件ug=ug/100;P=-m*ug;num=size(P, 1);c 二36612; % 阻尼k 1=2350000; %初始刚度k2二kl*0・2;a=zeros(l t num) ;v=zeros(l, num) ;x=zeros(l, num);%加速度速度位移Ei=zeros(l, num);Ek二zeros(1, num);Ed二zeros(1, num); Eh二zeros (1, num)输入能动能阻尼耗能滞回耗能EI=zeros(l t num); EK二zeros (1, num) ;ED=zeros(l, num) ;EH=zeros(l, num) ;%累积的各种能量t imc=zeros(l, num);a(l)=P (l)/m;hfl=zeros (1, num);h=0. 02;%地震波采样间隔Xy=0.01;%屈服位移pxmax=0;nxmax-0;pd=0;%双线型滞冋模型折线的标识，0农示弹性，1表示正向弹犁性，2表示反向弹性，-1 表示反向弹塑性，-2表示正向弹性for ii=l:numif pd==0 %弹性阶段k二kl;if x(ii)>Xypd=l;b=[(a(ii)-a(ii-l))/6/h a(ii~l)/2 v(ii-l) x(iiT)-Xy];% 拐点处理d=roots(b);for j=l:length(d)if isreal (d(j))==ldth=d(j);endendhp二h-dth;vp=v(ii-l)+a(ii-l)*dth+((a(ii)-a(ii-l))*dth*2/h/2);ap=a(i i-l) + (a (i i)-a(ii-l)) *dth/h;pp=m*ap+c*vp+kl*Xy;kd=k2+3*c/hp+6*m/hp/hp;dtpd=P(i i +1) -pp+m* (6*vp/hp+3*ap) +c* (3*vp^-hp*ap/2); dtx=dtpd/kd;dtv=3*dtx/hp-3*vp-hp*ap/2;x(ii)=Xy+dtx;v(ii)=vp^dtv;dtf=k2*dtx;hfl(ii)=kl*Xy+dtf;a(ii) = (P(ii)-hfl (i i+1)-c*v(i i))/m;elseif x(ii)<-Xypd=-l;b=[(a(ii)-a(ii-l))/6/h a(ii-l)/2 v(iiT) x(iiT)+Xy];% 拐点处理d=roots(b);for j=l:length(d)if isreal (d(j))==ldth=d(j);endendhp=h-dth;vp=v(i i 1)・a(i i 1)*dth+ ((a(ii)-a(ii-l))*dth"2/h/2); ap a(ii-l) + (a(ii)-a(ii-l))*dth/h;pp=m*ap+c*vp_k1*Xy; kd=k2+3*c/hp+6*m/hp/hp;dtpd=P (i i + 1)-pp+m*(6*vp/hp+3*ap)(3*vp^hp*ap/2);dtx=dtpd/kd;dtv-3*dtx/hp-3*vp-hp*ap/2;x(ii)=-Xy+dtx： v(i i)=vp^dtv;dtf=k2*dtx; hfl(ii)=-kl*Xy+dtf;a(i i) = (P(i i)-hf 1 (i i)-c*v(ii))/m;endendif pd==l %正向弹塑性k=k2;if v(ii)<0pd=2;b=[(a(ii)-a(ii-l))/2/h a(ii-l) v(ii-l)];% 拐点处理d二rools(b);for j=l:length(d)if isreal (d(j))==ldth=d (j);endendhp=h-dth;xp=x(iiT)+v(i i-1) *dth+a (i i-1) *dth*2/2+ (a (i i) -a (i i-l))*dth*3/h/6; ap=a (i i-1)+ (a (i i) -a (i i-1)) *dth/h;pp-m*ap+c*O+k1*Xy+k2* (xp-Xy);kd=kl+3*c/hp+6*m/hp/hp:dtpd=P(ii+l)-pp+m*(6*0/hp+3*ap)+c*(3*0+hp*ap/2);dtx=dtpd/kd;dtv=3*dtx/hp-3*0-hp*ap/2;x(i i)=xp+dtx;v(ii)=O+dtv;dtf=kl*dtx;hfl(ii)=kl*Xy+k2*(xp-Xy)+dtf;a(ii) = (P(ii)-hfl (ii)-c*v(ii))/m;pxmax=xp;endif pd==2 %反向弹性k=kl;if x (i i)>pxmaxpd=l; ■b=[(a(ii)-a(ii-l))/6/h a(ii-l)/2 v(ii-l) x(ii-l)-pxmax]:% 拐点处理d=roots(b);for j=l:length(d)if isreal (d(j))==I dth=d(j);endendhp=h-dth; vp=v(ii-l)+a(ii-l)*dth+((a(ii)-a(ii-l))*dth"2/h/2); ap=a(ii-l)+(a(ii)-a(ii-l))*dth/h;pp-m*ap+c*vp+k1*Xy+k2*(pxmax-Xy); kd=k2+3*c/hp+6*m/hp/hp;dtpd=P(ii+1)-pp+m*(6*vp/hp+3*ap)+c*(3*vp+hp*ap/2); dtx=dtpd/kd;d t v=3*d1x/hp-3*vp-hp*ap/2;x(i i)二pxmax+dtx;v(ii)=vp+dtv;dtf=k2*dtx;hfl (i i)=kl*Xy+k2*(pxmax-Xy) +dtf;a(ii) = (P(ii)-hfl (ii)-c*v(ii))/m;el seif x (i i)<(pxmax~2*Xy)pd=-l;―.・pu og 、(( =)A *3—(二二Jq —(二)d )H(二)P,±P +AX盖〒Ax*上—殳艾舌 X F (二二jqEp*znp 2p+dA:=(=)A-xlp+AX*zlxpulxdu(二)x,z 、dB*dqld 家gldv、xlp^HA4p&<p d "P H X 4p-(z 、dp*dfdA*g)芯+(du*g+dq、dA*9)*u+dd —(7二)dupdpp -dq'dq'lu 关 9+dqb苗+0上HP上-(xx*^IAX*>l —z>l*xeulxd) +d A *o +d e *u "dd2>-p*(〒二)E —(=)B)+U丄一)fde二z 'q 'Fe p 決((TT OF(U)P ))+±P*II )P +UI )AH D A召p —qudqp u a)pUO」(「)p*pTH((「)P 二p冒岛(P)£8U 上二丄joj-(q)s_oof p三P+O A=)Awp+dxH(二)x-z 、dydq —0^—dq、w'pdup,(z 、d¥dq+o為)*¥(dp*e+dq/0*9)*lu+dd —(I +二)dupdlp‘ d<dq 、ul*9+d<。