MATLAB与工程应用 第7章 动力学与振动

subplot(2,1,1); plot(t,x(:,1),lintype(2*(i-1)+1:2*i)); hold on subplot(2,1,2); plot(x(:,1),x(:,2),lintype(2*(i-1)+1:2*i)); hold on end subplot(2,1,1);
2 2
dτ dx 2 4π 2 2 = x1 x 4 − 2 dτ x1
2
dx3 = x4 dτ dx 4 2 x2 x4 =− dτ x1
d r 4π 2 d 2a dr da da − r = − 2 KK r 2 + 2 =0 2 dτ dτ dτ r dτ dτ
带入前面的微分方程组，可得四个一阶微分方程。 带入前面的微分方程组，可得四个一阶微分方程。
Logspace, rad2deg, loglog, semilogx see help
编制执行文件frequency72.m，求频率响应。 一.编制执行文件 编制执行文件 ，求频率响应。
m = 1; zeta = 0.1:0.1:1; k = 1; wn = sqrt(k/m); 10 w = logspace(-1,1,400); rad2deg = 180/pi; s = j*w; for cnt = 1:length(zeta) xfer(cnt,:)=(1/m) ./ (s.^2 + 2*zeta(cnt)*wn*s + wn^2); mag(cnt,:) = abs(xfer(cnt,:)); phs(cnt,:) = angle(xfer(cnt,:))*rad2deg; end for cnt = 1:length(zeta) figure(1) loglog(w,mag(cnt,:),'k-')
0
function xdot=ForceOcillation(t,x,dummy,zeta,Omega,x0)
xdot=[x(2);-2.0*zeta*x(2)-x(1)+x0*cos(Omega*t)];
为了获得频谱图,建立函数文件 为了获得频谱图 建立函数文件 AmplitudeSpectrum.m
第七章 动力学与振动
7.1轨迹 轨迹 7.2单自由度系统 单自由度系统 7.3多自由度系统 多自由度系统
7.1轨迹 轨迹 举例说明:重力场中有两个物体 其中质量为m2 重力场中有两个物体,其中质量为 举例说明 重力场中有两个物体 其中质量为 的物体固定,而质量为 的物体绕m2做平面圆 而质量为m1的物体绕 的物体固定 而质量为 的物体绕 做平面圆 周运动.做圆周运动的 做圆周运动的m1物体的轨道半径用变 周运动 做圆周运动的 物体的轨道半径用变 表示,角度用变量 表示. 量r表示 角度用变量 表示 表示 角度用变量a表示 m1 r a m2 两物体系统
由初始条件（位移和速度均为 时 建立执行文件menu72.m 由初始条件（位移和速度均为1时）建立执行文件
zeta=[0.1 1.0 5.0]; tspan=linspace(0,40,400); lintype=['-b' '--r' '---r']; for i=1:3
[t,x]=ode45('FreeOcillation',tspan,[1 1],[],zeta(i));
function[f,amplitude]=AmplitudeSpectrum(yy,Fs,Nstart,N); f=(Fs*(0:N-1)/N)*2.0*pi; amplitude=abs(fft(yy(Nstart:Nstart+N),N))/N;
采样速率30/6000=0.005,则采样频率 则采样频率 采样速率 1/0.005=200,这个频率远远超出了必须达到的 这个频率远远超出了必须达到的 采样频率,结果显示截短频谱图 结果显示截短频谱图,需设置 采样频率,结果显示截短频谱图,需设置 Nstart=3200,N=2^11=2048。 。 fft的应用见 的应用见Help 的应用见 编制执行文件menu72f.m 编制执行文件
程 序 运 行 结 果
7.2.4 线性系统的频率响应、阶跃响应及 线性系统的频率响应、 脉冲响应
单自由度振动系统的强迫振动微分方程可为： 单自由度振动系统的强迫振动微分方程可为：
通过LAPLACE变换，得传递函数： 变换，得传递函数： 通过 变换
Байду номын сангаас
其中： 其中： ω = k m n
c ζ = 2 mω n
mω n
阻尼因子 引入新变量转化状态空间方程形式： 引入新变量转化状态空间方程形式：
x1 = x dx x2 = dτ
dx1 = x2 dτ dx 2 3 = −2ζx 2 − x1 − αx1 + X 0 f (τ ) dτ
ζ =
c 2 mω n
7.2.2 线性系统的自由振动
一.运动微分方程 运动微分方程 当 α = F (τ ) = 0 时，得到线性振动系统的自由 振动方程。 振动方程。 d2x dx + 2ζ +x=0 2 dτ dτ 二.MATLAB求解 求解 对应的函数文件FreeOcillation.m 对应的函数文件
yy=x(:,1); N=2048;Nstart=3200;Fs=200; [f,Amplitude]=AmplitudeSpectrum(yy,Fs,Nstart,N); semilogy(f(1:40),2*Amplitude(1:40)); xlabel('Frequency'); ylabel('Amplitude'); title('Response spectrum of a linear system'); hold on subplot(2,1,1); xlabel('Time( \tau)'); ylabel('Displacement x( \tau)'); title('Response of a linear system'); hold on
title('SDOF frequency response magnitudes for zeta = 0.2 to 1.0 in steps of 0.2') xlabel('Frequency(rad/sec)') ylabel('Magnitude') grid hold on end hold off for cnt = 1:length(zeta) figure(2) semilogx(w,phs(cnt,:),'k-') title('SDOF frequency response phases for zeta = 0.2 to 1.0 in steps of 0.2') xlabel('Frequency(rad/sec)') ylabel('phase') grid hold on end hold off
程 序 运 行 结 果
7.2单自由度系统 单自由度系统 7.2.1概述 概述
一.力学模型 力学模型 X（t） （）
c K，a
F（t）=X（0）kf（t） （） （ ） （）
m
弹簧—质量 阻尼系统 弹簧 质量—阻尼系统 质量 其中：振体质量为m，弹簧的线性系数为k， 其中：振体质量为 ，弹簧的线性系数为 ，非线 性系数为a，阻尼系数为c，外力F（ ）。 性系数为 ，阻尼系数为 ，外力 （t）。
卫星绕地球转动时， 等于地球的质量 等于地球的质量， 卫星绕地球转动时，m2等于地球的质量， m1等于卫星的质量，r为卫星球心与地球 等于卫星的质量， 等于卫星的质量 为卫星球心与地球 球心间的距离。 球心间的距离。其运动轨迹由下列方程组 决定： 决定：
式中： 其中t是时间变量 是时间变量， 为 式中：τ = t / p ，其中 是时间变量，p为 物体在地球表面做圆周运动的周期。在地 物体在地球表面做圆周运动 球表面，r=6.373x106 m。
建立函数文件orbit.m 建立函数文件 function xd=orbit(t,x) xd=[x(2);x(1)*x(4)^2-4.0*pi^2/x(1)^2; x(4);-2.0*x(2)*x(4)/x(1)]; 三组初始条件(t=0)： 三组初始条件 ：
组 1 2 3
X1初始 X2初始 X3初始 X4初始 轨迹类型
zeta=0.4;Omega=3.0;x0=50; tspan=linspace(0,30,6000); options=odeset('RelTol',1e-8,'AbsTol',1e-8); lintype=['-b']; [t,x]=ode45('ForceOcillation',tspan,[0 0],options,zeta,Omega,x0); subplot(2,1,1); plot(t,x(:,1)); axis([0 30 -8 8]); hold on subplot(2,1,2);
function xdot=FreeOcillation(t,x,dummy,zeta)
xdot=[x(2);-2.0*zeta*x(2)-x(1)]; 三种阻尼系数（ ）阻尼系数为0.1时是欠阻尼 三种阻尼系数（1）阻尼系数为 时是欠阻尼 情况（ ）阻尼系数为1时是临界阻尼情况 情况（2）阻尼系数为 时是临界阻尼情况 （3）阻尼系数为 时是过阻尼情况 ）阻尼系数为5时是过阻尼情况
xlabel('Time( \tau)'); ylabel('Displacement x( \tau)'); title('Displacement as a function of( \tau)'); axis([0 40 -2.0 2.0]); text(2.7,-1.3,'阻尼系数=0.1'); text(3.6,-0.1,'1.0'); text(3.6,1.0,'5.0'); subplot(2,1,2); xlabel('Displacement'); ylabel('Velocity'); title('Phase portrait'); axis([-2.0 2.0 -2.0 2.0]); text(0.7,-1.25,'阻尼系数=0.1'); text(0.8,-0.65,'1.0'); text(0.8,0.1,'5.0');