matlab仿真自由落体

通信仿真技术与实践上机作业一

（实例 1.1）试对空气中在重力作用下不同质量物体的下落过程进行建模和仿真。已知重力加速度g=9.8m/s^2，在初始时刻t0=0s时物体由静止开始坠落。考虑空气阻力的影响。

（1）建立数学模型

质量为m的物体在自由坠落过程中受到竖直向下的恒定重力和向上的空气阻力f 的作用，由牛顿第二定律，我们知道，重力mg，加速度a以及物体质量m之间的关系是：

mg-f=ma

f=k*(v^2) k=空气阻力系数，为一恒定值

a=g-k(v^2/m)

（2）数学模型的解析分析

v(t)=at

s(t)=21at^2

（3）根据数学模型建立计算机仿真模型（编程）

将方程转换为一种在自变量（时间）上的“递推”表达式

v(t+dt)=v(t)+dv=v(t)+adt

s(t+dt)=s(t)+ds=s(t)+v(t)dt

（4）执行仿真和结果分析

% 自由落体.m

% 模拟受到空气阻力的小球

g=9.8; % 重力加速度

a=g;

m=10;

k=0.5; % 空气阻力系数

v=0; % 设定初始速度条件

s=0; % 设定初始位移条件

t=0; % 设定起始时间

dt=0.1; % 设置计算步长

N=20; % 设置仿真递推次数. 仿真时间等于N与dt的乘积

for f=1:N

v=v+a*dt; % 计算新时刻的速度

a=a-k.*(v^2)./m; % 空气阻力f=k*(v^2)/m

s(f+1)=s(f)+v*dt; % 新位移

t(f+1)=t(f)+dt; % 时间更新

end

% 作图: 受空气阻力落体结果与自由落体结果对比

t=0:dt:N*dt;

subplot(1,2,2)

plot(t,s,'o');

xlabel('时间 t'); ylabel('位移 s');

legend('受空气阻力的落体');

运行得到的结果:

（5）仿真程序的功能扩展---以动态方式来观察物体坠落的过程

受到空气阻力落体动画.m

〔实例1.2〕对乒乓球的弹跳过程进行仿真。忽略空气对球的影响，乒乓球垂直下落，落点为光滑的水平面，乒乓球接触落点立即反弹。如果不考虑弹跳中的能量损耗，则反弹前后的瞬时速率不变，但方向相反。如果考虑撞击损耗，则反弹速率有所降低。我们希望通过仿真得出乒乓球位移随时间变化的关系曲线，并进行弹跳过程的“实时”动画显示。

（1）数学模型

首先对乒乓球弹跳过程进行一些理想化假设。设球是刚性的，质量为m，垂直下落。碰击面为水平光滑平面。在理想情况下碰击无能量损耗。如果考虑碰击面损耗，则碰击前后速度方向相反，大小按比例系0<=K<=1下降。在t时刻的速度设为v=v(t)，位移设为y=y(t)，并以碰击点为坐标原点，水平方向为坐标横轴建立直角坐标系。球体的速度以竖直向上方向为正方向。重力加速度为g=9.8m/s^2。

初始条件假设：设初始时刻t0=0球体的初始速度为v0=v(t0)，初始位移为y0=y(t0)。

受力分析：在空中时小球受重力F=mg作用，其中，g=-(dt/dv) 。则在t+dt时刻小球的速度为

v(t+dt)=v(t)-gdt 在t+dt时刻小球的位移为

y(t+dt)=y(t)+v(t)dt

在小球撞击水平面的瞬间，即y(t)=0的时刻，它的速度方向改变，大小按比例K衰减。

当K=1时，就是无损耗弹跳情况。因此，小球反弹瞬间（t+dt时刻）的速度为

v(t+dt)=-Kv(t)-gdt; 0<=K<=1

反弹瞬间的位移为

y(t+dt)=y(t)-Kv(t)dt=-Kv(t)dt

（2）仿真模型设计（程序）

从数学模型中可见，小球在空中自由运动时刻与撞击时刻的动力方程不同。通过小球所处位置（位移）是否为零可判定小球处于何种状态。程序文件代码如下。

乒乓球弹跳曲线.m

% 乒乓球弹跳曲线.m

g=9.8; % 重力加速度

v0=0; % 初始速度

y0=1; % 初始位置

m=1; % 小球质量

t0=0; % 起始时间

K=1; % 弹跳的损耗系数

N=5000; % 仿真的总步进数

dt=0.001; % 仿真步长

v=v0; % 初状态

y=y0;

for k=1:N

if (y>0)|(v>0) % 小球在空中的(含刚刚弹起瞬间)动力方程计算

v =v -g*dt;

y =y +v*dt;

else % 碰击瞬间的计算

y =y-K.*v*dt;

v =-K.*v-g*dt;

end

s(k)=y; % 将当前位移记录到s数组中以便作图

end

t=t0:dt:dt*(N-1);% 仿真时间长度

subplot(1,2,1);plot(t,s);

title('无损耗时的弹跳');

xlabel('时间 t');

ylabel('位移 y(t)');

axis([0 5 0 1.1]);

g=9.8; % 重力加速度

v0=0; % 初始速度

y0=1; % 初始位置

m=1; % 小球质量

t0=0; % 起始时间

K=0.85; % 弹跳的损耗系数

N=5000; % 仿真的总步进数

dt=0.001; % 仿真步长

v=v0; % 初状态

y=y0;

for k=1:N

if (y>0)|(v>0) % 小球在空中的(含刚刚弹起瞬间)动力方程计算 v =v -g*dt;

y =y +v*dt;

else % 碰击瞬间的计算

y =y-K.*v*dt;

v =-K.*v-g*dt;

end

s(k)=y; % 将当前位移记录到s数组中以便作图

end

t=t0:dt:dt*(N-1);% 仿真时间长度

subplot(1,2,2);

plot(t,s);

title('有损耗的弹跳曲线');

xlabel('时间 t');

ylabel('位移 y(t)');

axis([0 5 0 1.1]);