简介这份资源是一套基于MATLAB的猎鹰9号火箭发射与着陆仿真项目面向具备一定动力学与控制理论基础、希望借助编程实践理解航天工程的高年级本科生、研究生及工程爱好者。项目围绕火箭动力学建模展开涵盖重力、空气阻力与推力等因素并采用数值积分方法求解运动方程同时涉及发射阶段的推力优化与姿态控制以及着陆阶段的导航制导与状态反馈、PID控制器设计帮助读者把牛顿第二定律转化为可运行的仿真程序。压缩包共5个文件包含3个m脚本、1个docx说明文档和1个md文件整体约21KB脚本分别对应空气阻力计算、助推器着陆与发射仿真等模块文档则阐述模型假设、控制系统结构与仿真步骤。目前已有1152人学习下载适合作为控制系统课程设计或航天动力学入门的开源参考读者可据此运行、分析并修改代码加深对火箭飞行力学与MATLAB工程仿真的理解。1. 火箭升空 MATLAB 代码从一条“我设计、编写并求解”的仿真链路说起很多人第一次看到“火箭升空 MATLAB 代码”这个标题会以为它只是画一条漂亮的抛物线。真动手才发现难点根本不在画图而在把推力、质量变化、气动阻力、重力转弯这几件事同时塞进一个常微分方程组里还要让求解器不崩。我最初做这套仿真时目标很朴素让一枚两级火箭从发射台起飞一级分离、二级入轨再让一级返回着陆。听起来像游戏但每一步都对应真实的物理约束。这套代码适合两类人一类是想用 MATLAB 练手动力学仿真的学生或工程师另一类是已经会ode45但不知道如何把多级火箭、变质量系统和着陆段拼成完整链路的人。它不追求工程级精度但足够让你看清每个参数如何影响弹道也足够让你在调参时体会什么叫“差之毫厘谬以千里”。2. 把火箭飞行动力学拆成可求解的方程组状态量、力与质量2.1 为什么选二维纵向平面模型而不是六自由度真实火箭飞行是六自由度问题但入门仿真如果一上来就写六自由度大概率会在姿态动力学和坐标系变换里翻车。我一般建议先用二维纵向平面模型把核心矛盾跑通火箭只在发射平面内运动位置用高度和水平距离描述速度分解为垂直和水平分量。这样状态量只有六个水平位置、垂直位置、水平速度、垂直速度、剩余质量、飞行路径角。飞行路径角这个量很关键它决定了推力方向与速度方向的夹角直接影响重力转弯的效率。选二维模型的另一个理由是求解速度。ode45在六自由度下容易因为姿态角高频振荡而被迫缩小步长仿真一次要等很久。二维模型下同样的物理过程几秒就能跑完方便你快速迭代参数。代价是忽略了滚转、偏航和横向风但对于理解“火箭为什么能入轨”这件事纵向平面已经足够。2.2 变质量系统的推力与质量流量建模火箭和普通抛体最大的区别是质量在变。推力来自燃料燃烧后高速喷出的气体质量流量和推力通过比冲关联。我通常用下面这组关系% 变质量火箭纵向平面动力学模型 % 状态向量 y [x, h, vx, vh, m, gamma] % x: 水平距离(m), h: 高度(m) % vx: 水平速度(m/s), vh: 垂直速度(m/s) % m: 当前质量(kg), gamma: 飞行路径角(rad) function dydt rocketDynamics(t, y, params) x y(1); h y(2); vx y(3); vh y(4); m y(5); gamma y(6); % 当前速度大小 v sqrt(vx^2 vh^2); % 大气密度指数模型海平面密度1.225 kg/m^3标高8500 m rho 1.225 * exp(-h / 8500); % 气动阻力阻力系数Cd参考面积A D 0.5 * rho * v^2 * params.Cd * params.A; % 推力一级工作时按恒定推力燃尽后切换二级 if t params.t_stage1 T params.T1; mdot -params.T1 / (params.Isp1 * params.g0); elseif t params.t_stage2 T params.T2; mdot -params.T2 / (params.Isp2 * params.g0); else T 0; mdot 0; end % 推力方向沿飞行路径角简化处理 Tx T * cos(gamma); Th T * sin(gamma); % 阻力方向与速度相反 if v 0 Dx -D * vx / v; Dh -D * vh / v; else Dx 0; Dh 0; end % 重力 g params.g0 * (params.Re / (params.Re h))^2; % 加速度 ax (Tx Dx) / m; ah (Th Dh) / m - g; % 飞行路径角变化率由速度方向变化决定 dgamma (ah * vx - ax * vh) / (v^2 1e-6); dydt [vx; vh; ax; ah; mdot; dgamma]; end这段代码里mdot是质量流量负号表示质量减少。Isp1和Isp2分别是两级发动机的比冲g0是海平面重力加速度。推力方向我简化成沿飞行路径角这在重力转弯段是合理近似但在主动段早期会有偏差。如果你要做更精确的仿真可以把推力方向改成由俯仰程序角控制那需要额外引入姿态方程。参数方面Cd和A决定阻力大小。对于细长体火箭跨声速段Cd会急剧上升但入门仿真用常数Cd也能看出趋势。Re是地球半径用于计算重力随高度的衰减。这些参数没有唯一标准值你可以根据公开的运载火箭量级自己设定重点是理解每个参数对弹道的影响方向。2.3 用 ode45 求解并处理事件触发ode45是 MATLAB 里最常用的非刚性求解器适合这类动力学问题。但直接调用它会遇到两个问题一是火箭落地或入轨后求解器还在算二是级间分离时刻需要切换参数。我的做法是用事件函数检测高度过零和速度方向变化% 事件函数检测着陆和入轨条件 function [value, isterminal, direction] rocketEvents(t, y, params) h y(2); v sqrt(y(3)^2 y(4)^2); % 事件1高度降到0着陆或坠毁 value(1) h; isterminal(1) 1; direction(1) -1; % 只检测下降穿越 % 事件2速度达到入轨速度约7800 m/s value(2) v - 7800; isterminal(2) 1; direction(2) 1; % 只检测上升穿越 end事件函数返回的value在过零时触发isterminal为1表示触发后停止积分direction控制检测方向。这样你就不用担心求解器在火箭落地后继续算出一堆无意义的数据。调用方式options odeset(Events, (t,y) rocketEvents(t,y,params), ... RelTol, 1e-6, AbsTol, 1e-8); [t, y, te, ye, ie] ode45((t,y) rocketDynamics(t,y,params), ... [0, 1000], y0, options);RelTol和AbsTol是相对和绝对容差。火箭动力学里速度变化剧烈容差放太宽会导致能量不守恒弹道明显偏离。我一般用1e-6和1e-8再小会显著拖慢速度。ie返回触发的事件编号你可以据此判断是入轨还是着陆。3. 一级返回与着陆段从弹道仿真到推力反推3.1 着陆段为什么不能沿用上升段模型上升段火箭质量大、推力方向基本沿速度方向重力转弯自然发生。但一级返回时火箭已经分离质量小得多而且需要发动机反推减速推力方向与速度方向相反。如果直接沿用上升段代码你会发现火箭要么直接砸地要么速度发散。原因是着陆段需要精确控制推力大小和方向让垂直速度在触地前降到接近零。我通常把着陆段单独建模状态量简化为高度、垂直速度和质量。水平方向假设已经通过栅格翼调整到接近零只考虑垂直方向的一维运动。这样问题从六状态降到三状态求解更稳定也更容易调参。3.2 推力反推的两种实现开环程序与闭环控制开环程序就是预先设定推力随时间的变化曲线比如“最后10秒推力从额定值线性降到零”。这种做法简单但对初始条件敏感初始高度或速度稍有偏差就会着陆失败。闭环控制则根据当前高度和速度实时计算所需推力常用的是比例-微分控制% 着陆段闭环推力控制 function T landingThrust(h, vh, m, params) % 目标在高度h_target处垂直速度降到v_target h_target 0; v_target -2; % 允许轻微下降速度避免悬停 % 比例项高度偏差 Kp params.Kp; % 微分项速度偏差 Kd params.Kd; % 期望加速度让速度按线性规律趋近目标 a_des Kp * (h - h_target) Kd * (vh - v_target); % 所需推力克服重力并产生期望加速度 g params.g0; T m * (g a_des); % 推力限幅不能超过发动机最大推力也不能为负 T max(0, min(T, params.T_max)); endKp和Kd是控制增益需要根据火箭质量和发动机推力调。Kp太大会导致推力剧烈振荡Kd太大会让响应变慢。我一般先用小增益跑通再逐步加大到着陆精度满足要求。T_max是发动机最大推力限幅是为了防止计算出物理上不可能的推力值。闭环控制的优点是鲁棒性好初始条件在一定范围内变化都能着陆。缺点是增益调参需要经验而且如果高度或速度测量有噪声推力会抖动。入门仿真里假设测量无噪声所以闭环控制效果通常比开环好。3.3 用表格对比上升段与着陆段的关键参数参数上升段着陆段状态量6个位置、速度、质量、路径角3个高度、垂直速度、质量推力方向沿速度方向或程序角与速度方向相反质量变化燃料消耗质量减小燃料消耗质量减小主要外力推力、重力、气动阻力推力、重力求解器容差RelTol 1e-6RelTol 1e-8终止条件速度达到入轨速度高度降到零这张表不是让你照抄而是提醒你两段仿真的物理假设不同代码不能简单复用。我见过有人把上升段代码直接改个推力方向就用来做着陆结果求解器在触地前一步崩溃血泪经验。4. 避坑与排查火箭仿真里最容易翻车的五个地方4.1 现象求解器报“无法满足积分容差”原因通常是动力学方程在某个时刻出现奇点比如速度为零时计算阻力方向或者质量耗尽后mdot仍然非零导致质量为负。解决方法是给速度加一个小量1e-6避免除零同时在质量流量计算里加判断当质量小于干重时强制mdot 0。4.2 现象火箭飞出发射平面后高度突然变成负值这是坐标系定义混乱导致的。有人用高度作为状态量但在计算重力时用了地心距两者差一个地球半径。解决方法是统一用高度重力加速度按g0 * (Re/(Reh))^2计算不要混用。4.3 现象级间分离后速度出现跳变原因是分离时刻质量突变但速度没有相应调整。真实分离是爆炸螺栓断开质量瞬间减少速度不变。如果你在代码里直接改质量而不改速度加速度会突变求解器可能误判。正确做法是在事件触发时重新初始化状态向量保持速度连续只更新质量。4.4 现象着陆段推力振荡火箭反复上下这是控制增益过大导致的。Kp太大时高度偏差被放大成巨大的推力指令火箭减速过头又上升形成极限环。解决方法是先减小Kp增大Kd让速度变化更平滑。如果仍然振荡检查推力限幅是否合理T_max太小会导致推力饱和控制失效。4.5 现象仿真结果和公开弹道数据对不上先别怀疑代码检查单位。我见过有人把比冲的单位当成秒但公式里用了米每秒结果推力大了十倍。比冲的定义是Isp T / (mdot * g0)单位是秒。如果你用T Isp * mdot * g0g0必须用9.81。另一个常见问题是角度单位MATLAB三角函数默认弧度如果你输入的是度弹道会完全不对。5. 让仿真更可信能量校验、参数扫描与一个实用技巧仿真跑通只是第一步怎么知道结果可信我习惯做两件事能量校验和参数扫描。能量校验的思路是在无阻力、无推力段机械能应该守恒。你可以把阻力和推力设为零让火箭做抛体运动然后计算动能加势能随时间的变化。如果能量漂移超过1%说明求解器容差太宽或方程有误。这个检查能抓住大部分符号错误和单位错误。参数扫描则是把比冲、推力、阻力系数等参数在一定范围内变化观察入轨速度和着陆精度的敏感度。我通常用for循环批量跑把结果存成矩阵再用imagesc画敏感度热图。这样你能一眼看出哪个参数最值得优化。比如比冲提高5%对入轨速度的影响可能比推力提高10%还大。一个实用技巧是在ode45的选项里打开Stats它会输出求解过程中调用函数次数和步长信息。如果步长小得离谱说明方程刚性太强或存在不连续点。这时候可以考虑换ode15s或者把不连续点如级间分离拆成两段分别积分。我一般会在分离时刻强制重启求解器而不是让ode45自己跨过不连续点这样更稳。最后说个习惯每次改完参数先跑一遍无推力无阻力的抛体确认能量守恒再跑完整弹道。这个习惯帮我省了很多次“为什么弹道不对”的排查时间。火箭仿真没有后悔药但有一个可靠的校验流程能让你少走很多弯路。希望帮到你。本文还有配套的精品资源点击获取