简介一套面向航天动力领域工程师与Matlab仿真学习者的固体火箭发动机SRM模拟器源码包基于数值方法对点火、燃烧及推进过程进行建模与可视化分析。压缩包共含7个文件核心为Motor.m、Propellant.m、SRM.m三个Matlab脚本分别对应发动机总体、推进剂特性和主仿真流程另有sorbitol_coarse.br与sorbitol_fine.br两种燃烧剂数据文件以及README.md说明文档和.gitignore配置整体仅5KB结构精简适合快速上手。目前已有69人学习下载。资源覆盖几何模型建立、燃烧室与喷嘴分析、火药燃烧特性计算、运动方程求解及参数优化等典型环节通过阅读代码可理解固体火箭发动机仿真从数据输入到结果输出的完整实现路径。虽然包体小但提供了清晰的模块划分和基础数据便于初学者在Matlab中直接运行调试也可作为进一步开发复杂仿真模型与进行敏感性分析的起点对相关课程设计或预研项目具有实用参考价值。 固体火箭发动机的“性格”其实很直接——药柱点着之后几乎不受外部控制压强、推力从点火到熄火全凭设计参数说了算。所以搞固发的人比搞液体的人更依赖仿真液体发动机可以随时关阀门固体发动机点着就只能等它烧完。想在不花试车钱的情况下先把发动机的“脾气”摸清楚最常用的办法就是建一个Matlab里的固体火箭发动机模拟器。这个模拟器听起来高大上拆开看核心就是一个零维内弹道计算程序把推进剂燃烧产气、喷管排气、燃烧室自由容积这几个物理过程用常微分方程组描述然后用ode45这类积分器跑一遍就能得到燃烧室压强、推力、燃面退移随时间的变化。非常适合刚接触推进系统仿真的学生、做方案设计的工程师以及想给课程加一个硬核案例的高校老师。我在实作中把模型、代码、调参经验和踩坑记录完整走了一遍下面直接分享整个设计思路和Matlab实现细节。1. 为什么用Matlab给固体火箭发动机做仿真1.1 这个模拟器解决什么问题固体火箭发动机的研制流程里地面试车是绕不开的环节但试车成本高、周期长、测试项目多不可能每个设计想法都拉去点火。设计初期真正需要的是一个能快速回答“如果把喉径改小5%燃烧室压强会不会超限”“药柱肉厚增加20%工作时间大概会拉长多少”这类问题的工具。这个模拟器干的就是这件事把发动机从点火到熄火的压强、推力、燃面变化过程算出来让设计人员在投料加工之前先评估一轮方案可行性。除了工程预研教学场景也很实用。很多推进系统课程里学生对“平衡压强”“燃速压强指数”“喉部塞塞”的理解停留在公式推导层面一旦看到实际压强曲线在点火瞬间冲高、随后缓慢爬升的趋势概念一下子就串起来了。这个项目作为课程设计或毕业设计的入门题目尺度刚刚好数学上不复杂物理上覆盖了燃烧、热力学和气体动力学又能充分练习Matlab的数值求解和可视化能力。1.2 为什么选Matlab而不是其他工具有人可能会问同样的事用Python或C写不也行吗对但Matlab在原型验证阶段有三个不可替代的优势。第一ode系列求解器封装成熟从非刚性到刚性、从事件检测到参数调优都是现成的不用自己写数值积分算法。第二脚本式工作流迭代极快改一个燃速系数、重新跑一遍、立刻画图看对比整个过程可以控制在几十秒内这在方案扫参阶段非常舒服。第三工具箱配合度高后面如果想把模型接到Simulink里做控制系统仿真或者用App Designer做个参数输入面板基本是无缝衔接。Matlab的劣势是商业授权和运行效率但这个量级的模型计算量非常小性能根本不是瓶颈。对我个人而言在方案探索期“跑得快”比“跑得省”重要得多所以我最终选择了Matlab作为实现平台。1.3 模型层级怎么选固体火箭发动机内部流动实际上是一个三维、非定常、带燃烧和两相流动的复杂过程但如果一上来就做CFD光网格划分和湍流模型调试就得折腾几周。模拟器采用零维集总参数模型先抓住主要矛盾。模型层级基本假设计算成本适用阶段零维集总参数燃烧室内压强、温度处处均匀极低秒级完成方案论证、参数扫描、迭代设计一维准一流动沿轴向有压强和流速梯度中等分钟级完成初步性能评估、侵蚀燃烧估算二维/三维CFD空间分布完全解析高小时到天级完成关键部位详细分析、点火过程三维仿真这里选零维模型的原因很实际内弹道特性主要取决于燃烧室质量守恒和喷管排气能力空间分布对总压和总推力的影响在工程精度上可以后置修正。等到设计方案初步收敛再往一维或三维深化效率会高很多。2. 核心建模原理与关键方程2.1 质量守恒整个内弹道的地基固体火箭发动机燃烧室内做的事情可以浓缩成一句话推进剂燃烧不断产生燃气燃气通过喷管不断排出去。燃烧室压强的变化本质上是燃气生成量和喷管排出量不平衡的结果。燃气生成量和排出量相等时燃烧室处于平衡状态压强基本稳定生成量大于排出量时压强上升反之则下降。用质量守恒方程写出来就是[ \frac{d(\rho_g V_c)}{dt} \dot m_{gen} - \dot m_{nozzle} ]其中 (\rho_g) 是燃气密度(V_c) 是燃烧室自由容积(\dot m_{gen}) 是燃气生成率(\dot m_{nozzle}) 是喷管质量流率。如果假设燃气为理想气体压强、密度、温度满足 (\rho_g P_c/(R_g T_c))并且把燃烧室温度近似看作常数就可以推导出压强变化率的显式表达式[ \frac{dP_c}{dt} \frac{R_g T_c}{V_c}(\dot m_{gen} - \dot m_{nozzle}) - \frac{P_c}{V_c}\frac{dV_c}{dt} ]注意等式右边的第二项这是自由容积变化带来的修正项。随着药柱不断烧去燃烧室内孔扩大自由容积持续增加这一项对压强的长期漂移有明显影响不能随手丢掉。很多初版代码跑出来的压强曲线出现莫名其妙的缓慢下降或漂移往往就是漏了这一项。2.2 燃速模型与推进剂参数固体推进剂的燃速是内弹道计算的灵魂。用得最多的是Vieille燃速定律[ r_b a P_c^n ]其中 (r_b) 是线燃速单位是m/s(P_c) 是燃烧室压强系数 (a) 主要取决于推进剂配方和初温压强指数 (n) 决定了燃速对压强的敏感程度。复合推进剂比如AP/HTPB体系的 (n) 一般在0.2到0.5之间(n) 越接近0燃速越不依赖压强发动机的工作稳定性越好。设计上通常希望 (n) 不要太高否则压强小幅波动会被燃速变化放大形成不稳定循环。燃速系数 (a) 的典型量级在 (10^{-5}) 量级单位换算成m/s/Pa^n后具体值受氧化剂粒径、配方组分比例影响很大。初算阶段可以直接取一个合理值后续用试车数据反推修正。有了燃速再结合当前燃面面积 (A_b) 和推进剂密度 (\rho_p)就得到燃气生成率[ \dot m_{gen} \rho_p A_b r_b ]这里 (A_b) 的确定依赖装药几何构型。端面燃烧药柱的燃面是个常数从头到尾不变内孔燃烧药柱的燃面会随着肉厚烧去逐渐增大。所以模拟器里要把“当前燃面半径”作为状态变量推进而不是拍脑袋给一个固定值。2.3 喷管流量与特征速度喷管在正常工作状态下处于塞塞状态喉部流速达到当地声速流量只跟燃烧室压强、喉部面积和燃气性质有关与下游环境压力无关这就是火箭发动机能产生推力的基础。零维模型中喷管流量的简化形式是[ \dot m_{nozzle} \frac{P_c A_t}{c^*} ]其中 (A_t) 是喉部面积(c^) 是特征速度。特征速度由推进剂热力学性质决定可以查表、用热力计算软件估算或者直接用经验值。典型复合推进剂的 (c^) 在1500到1650 m/s之间本文算例取1550 m/s。这个式子把喷管排气能力压缩成了一个参数非常方便但代价是丢掉了燃气组分变化和两相流损失等细节。在零维内弹道阶段这个精度足够用。2.4 稳态压强估算一个趁手的公式在写代码之前可以先推导一个很有用的公式稳态平衡压强。当系统达到平衡时燃气生成量等于喷管排出量[ \rho_p A_b a P^n \frac{P A_t}{c^*} ]整理后得到[ P_{eq} \left( \frac{\rho_p A_b a c^*}{A_t} \right)^{\frac{1}{1-n}} ]这个式子看着简单实际非常有用。给定一个装药设计可以立刻估算出平衡压强判断是否在结构允许范围或试验台测量范围内。每次跑仿真之前先用这个公式估算一下如果数值仿真结果和公式差得太远说明代码大概率有bug。我用一个典型的小型试验发动机做算例内孔燃烧管状药柱长度0.5 m药柱外半径0.05 m初始内孔半径0.015 m推进剂密度1800 kg/m³燃速系数 (a3.5\times10^{-5})压强指数 (n0.35)喉部直径0.02 m特征速度1550 m/s。初始燃面面积约0.047 m²喉部面积约(3.14\times10^{-4}) m²。代入稳态公式[ P_{eq} \left( \frac{1800 \times 0.047 \times 3.5\times10^{-5} \times 1550}{3.14\times10^{-4}} \right)^{\frac{1}{0.65}} \approx 2.5 \text{ MPa} ]算出来大概是2.5 MPa这就是这台发动机的点火后平衡压强。整个计算只需要一分钟而且完全可以在写仿真代码之前完成作为检验数值结果的标尺。3. Matlab实现与核心代码3.1 参数结构体与主程序流程Matlab里组织参数最推荐用结构体所有物理参数集中放在一个变量里ODE函数和事件函数都能直接引用避免全局变量满天飞。主程序的任务是设置参数、定义求解器选项、调用ode45、画图。% 主程序固体火箭发动机零维内弹道模拟 clear; clc; % 定义发动机与推进剂参数 prm.rho_p 1800; % 推进剂密度kg/m3 prm.r_in 0.015; % 初始内孔半径m prm.R_out 0.05; % 药柱外半径m prm.L_grain 0.5; % 药柱长度m prm.e prm.R_out - prm.r_in; % 肉厚m prm.a 3.5e-5; % 燃速系数m/s/Pa^n prm.n 0.35; % 燃速压强指数 prm.A_t pi * 0.01^2; % 喉部面积m2 prm.c_star 1550; % 特征速度m/s prm.R_g 320; % 燃气气体常数J/(kg*K) prm.T_c 2800; % 燃气温度K prm.V0 0.002; % 初始自由容积m3 % 初始条件压强给一个小初值肉厚从0开始 y0 [1e5; 0]; % 求解器选项开启事件检测检测药柱烧完 opts odeset(Events, (t,y) burn_event(t, y, prm)); % 求解 0~20s实际会由事件函数提前终止 [t, y] ode45((t,y) motor_ode(t, y, prm), [0 20], y0, opts); % 从结果中提取压强、肉厚、燃面半径 Pc y(:,1); w y(:,2); r_burn prm.r_in w; % 估算推力曲线推力系数近似取1.5 CF 1.5; F CF * Pc * prm.A_t;这里有几个细节值得注意。初始压强给了1e5 Pa而不是0原因后面会专门说。求解时间跨度写的是0到20 s但事件函数会在药柱烧完时自动终止求解所以实际运行时间取决于肉厚和燃速大约几秒量级。3.2 核心ODE函数燃烧、排气和几何退移ODE函数是整个模拟器的心脏。状态向量只包含两个变量燃烧室压强 (P_c) 和已烧去肉厚 (w)。用肉厚作为状态变量而不是直接计算当前燃烧半径好处是数值积分天然保证几何退移的连续性不会因为压强振荡导致燃面半径跳变。function dydt motor_ode(t, y, prm) % y(1): 燃烧室压强Pa % y(2): 已烧去肉厚m Pc y(1); w y(2); % 如果药柱已烧完直接返回零变化 if w prm.e dydt [0; 0]; return; end % 当前燃面半径以及燃面面积内孔表面 两个端面 r_burn prm.r_in w; A_b 2 * pi * r_burn * prm.L_grain 2 * pi * r_burn^2; % 燃速与燃气生成率 r_b prm.a * Pc^prm.n; m_gen prm.rho_p * A_b * r_b; % 喷管流量 m_nozzle Pc * prm.A_t / prm.c_star; % 当前自由容积及自由容积变化率 V_c prm.V0 pi * (r_burn^2 - prm.r_in^2) * prm.L_grain; dVc_dt 2 * pi * r_burn * prm.L_grain * r_b; % 压强变化率包含自由容积变化修正 dPc_dt (prm.R_g * prm.T_c / V_c) * (m_gen - m_nozzle) ... - (Pc / V_c) * dVc_dt; % 肉厚增长率就是燃速 dw_dt r_b; dydt [dPc_dt; dw_dt]; end代码逻辑简单但每一步都有明确的物理对应燃面面积驱动燃气生成压强驱动喷管流量两者的差值加上体积变化项共同决定压强变化率。这正好呼应了前面质量守恒方程的推导代码和物理方程是一一对应的。事件函数用来检测药柱烧完的时刻这样不需要手动估算积分终点。function [value, isterminal, direction] burn_event(t, y, prm) % 剩余肉厚从正变为0时触发终止 value prm.e - y(2); isterminal 1; direction 0; end3.3 结果可视化与曲线解读模拟跑完后用两个子图分别展示燃烧室压强和推力的时间历程。figure; subplot(2,1,1); plot(t, Pc/1e6, LineWidth, 1.5); xlabel(时间 (s)); ylabel(燃烧室压强 (MPa)); grid on; title(燃烧室压强-时间曲线); subplot(2,1,2); plot(t, F/1000, LineWidth, 1.5); xlabel(时间 (s)); ylabel(推力 (kN)); grid on; title(推力-时间曲线);看一下曲线形态一般来说会呈现三个典型阶段点火初期燃气生成率快速超过喷管排出率压强在极短时间内冲高并出现一个小幅过冲这就是“点火压强峰”随后进入平衡段压强缓慢爬升或下降具体趋势取决于燃面面积变化——内孔燃烧药柱的燃面随肉厚增大而增大所以平衡压强势头向上最后在药柱烧尽的临近时刻燃面开始快速萎缩燃气生成率下降压强迅速跌落回环境压力。这个“快升、缓变、速降”的形状几乎就是固体火箭发动机内弹道曲线的名片。看到自己代码生成这条曲线时整个物理过程的轮廓就一下子立起来了。4. 常见问题与调试技巧4.1 数值发散罪魁祸首往往是初始压强很多第一次写这个模型的人习惯把初始压强设成0结果ode45一跑就发散或者算出负压强。原因在于燃速定律 (r_b aP^n)当 (P0) 时燃速为0燃气生成率为0系统一开始就处在“没有燃气”的平静状态数值上很难被激励起来而且 (P^n) 在0附近梯度很大积分器容易失稳。解决方法是给一个合理的小初值比如1e5 Pa让系统前几步就能被“点燃”。更讲究的做法是先用稳态公式解出平衡压强然后直接用平衡压强的70%~80%作为初始值这样可以跳过点火建立阶段专注研究中段弹道特性。4.2 单位制混乱最容易踩的坑这个模型里牵涉的单位很杂压强有Pa和MPa混用、燃速系数有mm/s和m/s混用、特征速度有时从文献里拿到的单位不是国际单位制。我吃过一次亏因为把燃速系数从mm/s量级直接当成m/s代入算出来的平衡压强差了三个数量级折腾了半天才发现是单位问题。建议从头到尾统一用国际单位制压强Pa、长度m、时间s、质量kg、温度K。所有从文献或报告里拿到的参数在初始化阶段就换算好中间不要夹带非国际单位。如果习惯了用MPa看图那就在画图时再除以1e6不要在计算过程中混用。4.3 理论结果对不上仿真结果如果仿真跑出来的平衡压强和稳态公式对不上先不要怀疑公式按这个顺序排查先检查稳态公式里用的燃面面积是不是和ODE函数里一致再检查喷管面积有没有换算错误最后看自由容积修正项有没有漏掉。有一类隐蔽错误是燃面面积随肉厚的变化趋势写反了导致压强曲线方向完全反了。还有个常见问题是忽略自由容积变化项导致压强缓慢漂移。虽然这一项在数值上看起来不大但在长时间工作发动机里随着药柱逐渐烧空自由容积可能增加好几倍累积效应非常明显。我的建议是无论题目简单与否都把 (dV_c/dt) 这一项写全养成好习惯。4.4 调参心态与模型的定位整个模拟器最常见的三个问题、原因和解决办法可以汇总成下表问题现象可能原因解决办法求解器报错积分失败初始压强为0给1e5 Pa量级的初始值平衡压强与理论估算差很远单位制混用统一国际单位制压强曲线后期漂移明显漏掉自由容积变化修正项补上 (-P_c/V_c\cdot dV_c/dt)曲线出现高频振荡模型刚性增强换ode15s或ode23s求解烧完时间异常长燃速系数数量级出错核对a的单位和量级说到这里我这个模拟器已经能完成从参数输入到曲线输出的完整流程但它的定位始终是一个“快速评估工具”而不是精确预测工具。零维模型丢掉了侵蚀燃烧、热损失、两相流损失、点火瞬态等物理过程设计后期需要逐项往模型里补充这些修正。最后说一点我自己的体会。第一次跑通模拟器时我盯着压强曲线在点火瞬间冲上去又缓缓爬升那种感觉像是真的“看到”了发动机在试车台上的反应。这个项目的价值不只在于算出几条曲线而是把燃烧、流动、几何退移这些原本需要跨学科知识才能串起来的物理过程压缩成了一个可以反复试错的计算环境。后续扩展方向也很明确给模型加一维轴向流动项、估算侵蚀燃烧对燃速的增强效应、把结构传热和壳体温度场耦合进来或者用App Designer做一个带滑条和参数输入框的交互界面。固体发动机的设计迭代本质就是在小模型里先跑通趋势再去试车台验证细节、修正模型参数。仿真替代不了试车但能让你在试车之前心里先有一个底。本文还有配套的精品资源点击获取