计算火箭飞行空气阻力系数

计算火箭飞行空气阻力系数一问题描述及要求

小型火箭初始质量为1200千克，其中包括900千克燃料。火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生40000牛顿的恒定推力，当燃料用尽时引擎关闭。设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数记作k。

利用数据拟合,计算空气阻力系数k

1.观测数据保存于rocket.dat

2.图形显示拟合结果

二物理模型分析

本问题是变质量的受力问题，应当采用=dp

F

dt 而非

2

2

=

d x

F m

dt

计算。因此受力方程为：

2

dp dm dv

v m F mg kv dt dt dt

=+=−−

其中

015

m m t

=−，使用差分方法可以由速度求出对应的加速度，同时也由位移求出对应的速度，带入上述方程即可线型拟合出比例系数k。

三程序分析及设计

由上述受力分析可以得到下面的代码。先通过for循环构建新的差分处理后的数据，再使用拟合得出k值，并以此作图。

clc;

clear;

A=load('G:/rocket.dat');

for i=1:19

t(i)=(A(1,(i+1))+A(1,i))./2;

a(i)=(A(3,(i+1))-A(3,i))./(A(1,(i+1))-A(1,i));

velocitysquare(i)=((A(2,(i+1))-A(2,i))./(A(1,(i+1))-A(1,i))).^2;

F(i)=40000-9.8.*(1200-15.*t(i))-a(i).*(1200-

15.*t(i))+15.*sqrt(velocitysquare(i));

i=i+1;

end %差分法处理数据得到加速度以及相应时刻的速度

p=polyfit(velocitysquare,F,1);%拟合曲线

plot(velocitysquare,F,'*','color','r','linewidth',2);%绘制实际测量值

hold on;

plot(velocitysquare,p(1).*velocitysquare+p(2),'color','b','linewidth' ,1.2);%绘制拟合曲线

xlabel('v^2'); ylabel('F=40000-9.8×(1200-15×t)-a×(1200-15×t)+15×v'); text(0.5*max(100000),0.43*max(30000),sprintf('F=%5.4f v^2

disp(sprintf('k=%8.5f',p(1)));%输出k值

四程序运行结果及讨论

程序拟合得k= 0.33010。相应实测结果与以及由拟合结果反推的结果如下图1所示。

图1 v^2-f关系

解受力微分方程，将实际测量数据与k= 0.33010时的方程解析解绘制在同一张图上，如图2所示。相关微分方程代码见附录。

图2 拟合值与原始实际测量值

可以观察到，拟合值相较于实际测量值，无论是位移还是速度都相对偏大。这很有可能是由于出题人未意识到变质量受力问题不能简单使用dv F m dt

=公式计算。因此为了揣摩出题人的意图，使用2dv m F mg kv dt =−−重新计算，得到拟合k= 0.28980。其中v^2-f 关系如图3所示。

图3 v^2-f 关系

解受力微分方程，在k= 0.28980且使用2dv m F mg kv dt

=−−公式的情况下，得到图4拟合值与原始实际测量值的x-t 与v-t 关系图。

图4 拟合值与原始实际测量值

显然，图4相较于图2拟合更好。这说明了出题人存在考虑不周的可能。

综上所述，若考虑变质量，则k= 0.33010，否则k= 0.28980。现实情况下应更符合前者，但是此题数据更符合不考虑变质量对动量的影响，即后者。

附录 程序代码

1、2dv m

F mg kv dt

=−−情况下的拟合代码 clc;

clear;

A=load('G:/rocket.dat');

for i=1:19

t(i)=(A(1,(i+1))+A(1,i))./2;

a(i)=(A(3,(i+1))-A(3,i))./(A(1,(i+1))-A(1,i));

velocitysquare(i)=((A(2,(i+1))-A(2,i))./(A(1,(i+1))-A(1,i))).^2; F(i)=40000-9.8.*(1200-15.*t(i))-a(i).*(1200-15.*t(i)));

i=i+1; end %差分法处理数据得到加速度以及相应时刻的速度

p=polyfit(velocitysquare,F,1); %拟合曲线

plot(velocitysquare,F,'*','color','r','linewidth',2); %绘制实际测量值 hold on ;

plot(velocitysquare,p(1).*velocitysquare+p(2),'color','b','linewidth',1.2); %绘制拟合曲线

xlabel('v^2');

ylabel(' F=40000-9.8×(1200-15×t)-a×(1200-15×t) ');

text(0.5*max(100000),0.43*max(30000),sprintf('F=%5.4f v^2

+%5.4f',p));

disp(sprintf('k=%8.5f',p(1))); %输出k 值

2、2dp dm dv v m F mg kv dt dt dt

=+=−−并且k=0.33010情况下的受力微分方程求解 function dy=lalala(t,y)%微分方程函数

k=0.33010;

dy=zeros(2,1);%明确 dy 的维数，即一阶常微分方程的个数，且为列向量

dy(1)=y(2);

dy(2)=(-40000+9.8*1200-9.8*15*t+k*(y(2)^2)+15*y(2))/(15*t-1200);

clc;

clear;

tspan=[0,20];

[t,y]=ode45('lalala',tspan,[0 0]);