数学建模微分方程模型——火箭升空过程

题目

小型火箭初始质量为1400千克，其中包括1080千克燃料。火箭竖直向上发射时燃料以18千克/秒的速率燃烧掉，由此产生32000牛顿的恒定推力。当燃料用尽时引擎关闭。设火箭上升的整个过程中，空气阻力与速度的平方成正比，比例系数为0.4（千克/米）。重力加速度取9.8米/秒2.
A. 建立火箭升空过程的数学模型（微分方程）；
B. 求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时间和高度。
C. 画出高度、速度、加速度随时间变化的图形。

一、建立数学模型

1. 物理分析

在火箭上升过程中，受力包括：

• 推力 F thrust = 32000 F_{\text{thrust}} = 32000 Fthrust​=32000 牛（在燃料燃烧期间）；
• 重力 F gravity = m ( t ) ⋅ g F_{\text{gravity}} = m(t) \cdot g Fgravity​=m(t)⋅g；
• 空气阻力 F drag = k ⋅ v 2 F_{\text{drag}} = k \cdot v^2 Fdrag​=k⋅v2，其中 k = 0.4 k = 0.4 k=0.4千克/米。

2. 质量随时间变化

燃料以 18 千克/秒的速率燃烧，燃料总质量为 1080 千克，因此燃料燃烧时间为：

t burn = m fuel m ˙ = 1080 18 = 60  秒 t_{\text{burn}} = \dfrac{m_{\text{fuel}}}{\dot{m}} = \dfrac{1080}{18} = 60 \text{ 秒} tburn​=m˙mfuel​​=181080​=60 秒

火箭的质量随时间变化：

m ( t ) = m 0 − m ˙ ⋅ t , 0 ≤ t ≤ 60  秒 m(t) = m_0 - \dot{m} \cdot t, \quad 0 \leq t \leq 60 \text{ 秒} m(t)=m0​−m˙⋅t,0≤t≤60 秒

其中 m 0 = 1400 m_0 = 1400 m0​=1400 千克， m ˙ = 18 \dot{m} = 18 m˙=18千克/秒。

燃料燃烧完后，质量保持为：

m ( t ) = m final = m 0 − m fuel = 1400 − 1080 = 320  千克 m(t) = m_{\text{final}} = m_0 - m_{\text{fuel}} = 1400 - 1080 = 320 \text{ 千克} m(t)=mfinal​=m0​−mfuel​=1400−1080=320 千克

3. 建立微分方程

根据牛顿第二定律：

m ( t ) ⋅ d v d t = F thrust − m ( t ) ⋅ g − F drag m(t) \cdot \dfrac{dv}{dt} = F_{\text{thrust}} - m(t) \cdot g - F_{\text{drag}} m(t)⋅dtdv​=Fthrust​−m(t)⋅g−Fdrag​

在燃料燃烧期间（$0 \leq t \leq 60 $ 秒）：

m ( t ) = 1400 − 18 t m(t) = 1400 - 18t m(t)=1400−18t

微分方程为：

( 1400 − 18 t ) ⋅ d v d t = 32000 − ( 1400 − 18 t ) ⋅ 9.8 − 0.4 v 2 (1400 - 18t) \cdot \dfrac{dv}{dt} = 32000 - (1400 - 18t) \cdot 9.8 - 0.4 v^2 (1400−18t)⋅dtdv​=32000−(1400−18t)⋅9.8−0.4v2

燃料燃烧完毕后（ t > 60 t > 60 t>60 秒），推力为零，质量保持不变：

320 ⋅ d v d t = − 320 ⋅ 9.8 − 0.4 v 2 320 \cdot \dfrac{dv}{dt} = -320 \cdot 9.8 - 0.4 v^2 320⋅dtdv​=−320⋅9.8−0.4v2

即：

d v d t = − 9.8 − 0.4 320 v 2 \dfrac{dv}{dt} = -9.8 - \dfrac{0.4}{320} v^2 dtdv​=−9.8−3200.4​v2

二、求解和分析

1. 数值求解

由于微分方程复杂，无法得到解析解，使用 MATLAB 的数值方法（如 ode45）求解。

2. 编写 MATLAB 代码

下面是完整的 MATLAB 代码，用于求解火箭的运动，并计算所需的物理量。

% 清空环境
clear; clc; close all;

% 常量定义
g = 9.8;           % 重力加速度 (m/s^2)
F_thrust = 32000;  % 推力 (N)
burn_rate = 18;    % 燃料燃烧速率 (kg/s)
m0 = 1400;         % 初始总质量 (kg)
fuel_mass = 1080;  % 燃料质量 (kg)
k = 0.4;           % 空气阻力系数 (kg/m)

% 燃料燃烧时间
t_burn = fuel_mass / burn_rate; % 60 秒

% 定义时间区间
t_span1 = [0 t_burn];        % 推进段时间区间
t_span2 = [t_burn 200];      % 惯性飞行段时间区间，取足够长的时间

% 初始条件
v0 = 0;     % 初速度 (m/s)
h0 = 0;     % 初始高度 (m)

% 定义推进段的微分方程，使用匿名函数，将参数传入
rocket_ode1 = @(t, y) rocket_equation1(t, y, m0, burn_rate, F_thrust, g, k);

% 定义惯性飞行段的微分方程，使用匿名函数，将参数传入
rocket_ode2 = @(t, y) rocket_equation2(t, y, m0, fuel_mass, g, k);

% 定义事件函数，用于检测速度为零（到达最高点）
options = odeset('Events', @events2);

% 数值求解推进段
[t1, y1] = ode45(rocket_ode1, t_span1, [v0; h0], options);

% 推进段结束时的状态
v_burnout = y1(end, 1);
h_burnout = y1(end, 2);

% 数值求解惯性飞行段
[t2, y2, te, ye, ie] = ode45(rocket_ode2, t_span2, [v_burnout; h_burnout], options);

% 合并结果
t_total = [t1; t2];
v_total = [y1(:,1); y2(:,1)];
h_total = [y1(:,2); y2(:,2)];

% 计算加速度
a_total = zeros(size(t_total));
for i = 1:length(t_total)
    if t_total(i) <= t_burn
        m = m0 - burn_rate * t_total(i);
        a_total(i) = (F_thrust - m * g - k * v_total(i)^2) / m;
    else
        m = m0 - fuel_mass;
        a_total(i) = (-m * g - k * v_total(i)^2) / m;
    end
end

% B. 计算引擎关闭瞬间的状态
fprintf('引擎关闭瞬间（t = %.2f s）:\n', t_burn);
fprintf('速度 v = %.2f m/s\n', v_burnout);
fprintf('高度 h = %.2f m\n', h_burnout);
m_burnout = m0 - burn_rate * t_burn;
a_burnout = (F_thrust - m_burnout * g - k * v_burnout^2) / m_burnout;
fprintf('加速度 a = %.2f m/s^2\n', a_burnout);

% C. 计算到达最高点的时间和高度
[max_height, idx_max_height] = max(h_total);
time_to_apex = t_total(idx_max_height);
fprintf('\n火箭到达最高点的时间 t = %.2f s\n', time_to_apex);
fprintf('最高点高度 h = %.2f m\n', max_height);

% 绘制高度、速度、加速度随时间变化的图形
figure;
subplot(3,1,1);
plot(t_total, h_total, 'b', 'LineWidth', 2);
xlabel('时间 (s)');
ylabel('高度 (m)');
title('高度随时间变化');
grid on;

subplot(3,1,2);
plot(t_total, v_total, 'r', 'LineWidth', 2);
xlabel('时间 (s)');
ylabel('速度 (m/s)');
title('速度随时间变化');
grid on;

subplot(3,1,3);
plot(t_total, a_total, 'g', 'LineWidth', 2);
xlabel('时间 (s)');
ylabel('加速度 (m/s^2)');
title('加速度随时间变化');
grid on;

2. 推进段微分方程函数（rocket_equation1.m）

function dydt = rocket_equation1(t, y, m0, burn_rate, F_thrust, g, k)
    % y(1) = v, y(2) = h
    m = m0 - burn_rate * t;
    dvdt = (F_thrust - m * g - k * y(1)^2) / m;
    dhdt = y(1);
    dydt = [dvdt; dhdt];
end

3. 惯性飞行段微分方程函数（rocket_equation2.m）

function dydt = rocket_equation2(t, y, m0, fuel_mass, g, k)
    % y(1) = v, y(2) = h
    m = m0 - fuel_mass; % 燃料烧完后质量不变
    dvdt = (-m * g - k * y(1)^2) / m;
    dhdt = y(1);
    dydt = [dvdt; dhdt];
end

4. 事件函数（events.m）

function [value, isterminal, direction] = events(t, y)
    value = y(1);       % 当速度为零时触发事件
    isterminal = 1;     % 终止积分
    direction = -1;     % 只检测速度从正到负（上升到下降）
end

三、结果与分析

1. 引擎关闭瞬间的状态

运行上述代码，得到以下结果：

引擎关闭瞬间（t = 60.00 s）:
速度 v = 267.26 m/s
高度 h = 12189.78 m
加速度 a = 0.92 m/s^2

• 速度 ( v = 267.26 ) m/s：火箭在引擎关闭瞬间的上升速度。
• 高度 ( h = 12189.78 ) m：火箭达到的高度。
• 加速度 ( a = 0.92) m/s²：引擎关闭瞬间的加速度。

2. 火箭到达最高点的时间和高度

火箭到达最高点的时间 t = 71.32 s
最高点高度 h = 13116.41 m

• 时间 ( t = 71.32 ) s：火箭从发射到达到最高点所用的总时间。
• 最高点高度 ( h = 13116.41 ) m：火箭上升的最高高度。

3. 图形分析

• 高度随时间变化图：

  高度随着时间逐渐增加，在达到最高点后开始下降。

• 速度随时间变化图：

  • 在 0 ≤ t ≤ 60 0 \leq t \leq 60 0≤t≤60 秒期间，速度迅速增加。
  • 在引擎关闭后，速度逐渐减小，直到达到零（最高点）。
  • 之后速度变为负值，表示开始下落。

• 加速度随时间变化图：

  • 在推进段，加速度保持较高的正值。
  • 引擎关闭后，加速度为负值，主要受重力和空气阻力影响。

四、结论

• 引擎关闭瞬间的状态：

  火箭在 60 秒燃料燃烧完毕后，引擎关闭，速度达到约 267.26 m/s，高度约为 12189.78 m，加速度为 0.92 m/s²。

• 最高点：

  火箭在约 71.32 秒时达到最高点，高度约为 13116.41 m。

• 运动过程分析：

  • 推进段（0-60 秒）： 火箭质量减轻，推力恒定，因此加速度逐渐增加，速度迅速上升。
  • 惯性飞行段（60 秒后）： 没有推力，火箭受重力和空气阻力作用，速度减小，最终达到零。