简介面向固体火箭发动机设计与仿真领域的Matlab模拟器压缩包适合航天、机械及计算机仿真方向的研究者与工程师用于在物理试验前完成点火、燃烧、推进等环节的虚拟验证。包内共7个文件包含3个Matlab源程序覆盖主流程、发动机结构与推进剂模块2个br格式推进剂数据文件提供不同颗粒度燃烧参数另有说明文档与Git配置文件整体仅5KB结构精简。已有69人学习可用于快速搭建发动机几何模型、分析燃烧室与喷嘴流场、建立火药燃烧模型并求解运动方程涵盖从点火到推进的全过程仿真。借助Matlab可视化能力可直观观察温度、压力和流场变化支持参数优化与敏感性分析为设计验证和故障分析提供支撑对缩短研发周期、提升设计可靠性具有实用价值。 做固体火箭发动机内弹道仿真这件事我一开始并没有打算自己造轮子。当时手头有个小型固体火箭发动机的预研需求需要快速得到一组燃烧室压强曲线和推力曲线用来评估装药方案到底靠不靠谱。翻了一圈现成工具后我发现一个尴尬局面通用内弹道程序要么药型库固定想改几何参数得动源代码要么干脆是个黑盒算完给你一张曲线图中间压强怎么建立、燃面怎么发展、自由容积怎么变化全都看不见。用CFD又太重网格、湍流模型、两相流一轮下来项目周期根本扛不住。最后我决定在Matlab里手写一个固体火箭发动机模拟器基于零维内弹道模型把燃速模型、药柱几何、喷管流率和ODE求解完整串起来。跑通后十来分钟就能出一组结果而且每一步物理过程都能拆开看后来这个模拟器陪我改了好几轮方案也提前筛掉过两个参数上明显不合理的药型。这篇内容适合对固体火箭发动机有基本概念、想真正把零维内弹道代码跑通的工程师或学生我会把物理模型、代码实现、数值坑和验证方法一次讲清楚。1. 为什么自己写一个Matlab模拟器不直接套现成软件1.1 这个模拟器能算什么固体火箭发动机内弹道仿真的核心产出其实就是两条曲线燃烧室压强随时间的变化以及推力随时间的变化。这两条曲线几乎决定了发动机的所有主要性能——总冲、平均推力、最大压强、工作时间、装药是否合理。再往下还能导出燃面面积随时间的退化、推进剂剩余质量、点火瞬间的压强冲击高度。对预研阶段来说这个量级的信息已经完全够用。你不需要知道燃烧室内部哪里的流场有回流、哪里的温度更高因为零维模型本身就是把整个燃烧室当作一个充分混合的控制体压强在任意瞬间处处相等。这个假设听起来很粗暴但对绝大多数固体火箭发动机的设计迭代来说它给出的结果已经相当能打。用这套模拟器你可以做的事情包括比较不同药柱内径下的压强平台段评估喉部面积变大后对工作压强的影响估算点火瞬态的峰值压强计算总冲和平均比冲。这些都是在方案阶段必须回答的问题而且这些问题用手算只有稳态解看不到动态过程用CFD又像是在用大炮打蚊子。1.2 为什么选Matlab而不是专用软件有人会问NASA的CEA、商业软件里的固体发动机模块都能算为什么要自己在Matlab里折腾。我当时的判断很简单对比项Matlab模拟器专用内弹道程序通用CFD学习成本低核心方程能自己推导中需要熟悉输入格式高需要网格和模型经验可修改性完全开放改一行代码就换一种药型受程序架构限制改几何就要重新建模计算速度秒级到分钟级快小时级到天级结果透明度每个中间变量都能看只有最终曲线后处理复杂前期投入只花时间可能涉及授权人力物力都大Matlab还有个天然优势ODE求解器非常成熟数值刚性检测、事件触发、误差控制都有完整方案而且绘图、数据处理、参数扫描都在同一个环境里完成不用在Python、Origin、Excel之间来回倒腾。如果你在高校或研究所Matlab基本是标配不存在环境门槛。我把这套模拟器写成纯脚本加函数不依赖Simulink装个基础版Matlab就能跑。2. 零维内弹道模型的物理方程从质量守恒到压强微分2.1 核心闭环三个变量决定了燃烧室压强零维内弹道模型的根基就是燃烧室内的质量守恒。压强随时间的变化本质上由三个量在博弈推进剂燃烧产生的气体、喷管排出的气体、以及燃烧室自由容积的变化。写成微分方程就是dp_c/dt (R_gas * Tc / Vc) * (ṁ_gen - ṁ_nozzle) - (pc / Vc) * dVc/dt这里每一项都要拆开看。第一项里的R_gas是燃气的气体常数Tc是燃烧温度Vc是当前燃烧室自由容积。ṁ_gen是燃烧生成的气体质量流率正比于推进剂密度、当前燃面面积和燃速ṁ_nozzle是喷管排出的质量流率。最后一项里的dVc/dt是自由容积的变化率因为药柱不断烧掉内孔扩大留给气体的空间在不断变大。这个方程要说直观也直观生成比排出多压强升高自由容积增大压强倾向下降。两者达到平衡时dp_c/dt为零发动机就进入稳态工作段。这里我强调一点很多简化模型会忽略dVc/dt这一项但点火段和薄药柱的末段这一项的影响会被放大如果完全不考虑压强曲线会显得过于平直和实测对不上。2.2 燃速模型与药柱几何参数全在几何里燃速采用固体推进剂最经典的Vieille经验公式r_b a * pc^na是燃速系数n是压强指数。压强指数n的大小直接影响发动机的稳定性一般复合推进剂在0.3到0.5之间。指数越高压强波动越容易被放大所以设计时通常想办法压低n。参数单位是最大的坑。很多资料里给出a的单位是mm/s·MPa^-n但SI体系下压强单位是Pa直接代进去结果差好几个数量级。换算方法很简单如果查到的a_cmps是以mm/s和MPa为基准的转成SI要乘上(1e6)^(-n)因为1 MPa等于1e6 Pa。比如a1.2 mm/s·MPa^-0.4转成SI就是1.2e-3再乘以1e6的负0.4次方。这个换算我见过太多人栽跟头算出来的压强要么离谱地高要么离谱地低。药柱几何我以最常用的圆柱内孔药柱为例两端做阻燃包覆只有内孔表面参与燃烧。当前燃面直径d等于初始内径d_i加上两倍的已烧蚀厚度e燃面面积Ab π * d * L_grain。随着燃烧推进d变大Ab也跟着变大所以燃面是增面燃烧平衡压强会缓慢上升直到外层烧完。自由容积同步更新Vc V0 π/4 * (d^2 - d_i^2) * L_grainV0是装药前就存在的初始空腔容积也包括点火器空间。这个V0对点火瞬态影响极大后面我会专门说。2.3 喷管流率与特征速度c*喷管排出项写的是ṁ_nozzle pc * A_t / c*其中A_t是喉部面积c*是特征速度。这里没有真的去解喷管内流动而是用c*这个参数把燃烧室到喉部之间的能量转化打包描述。c*可以由燃烧温度、燃气常数、比热比算出c* sqrt(R_gas * Tc / γ) / sqrt((2/(γ1))^((γ1)/(γ-1)))工程上更常见的是直接用推进剂手册里的实验值比如典型的AP/HTPB复合推进剂c*大约在1500到1650 m/s。c*只反映燃烧气体的能量水平跟喷管扩张比和出口条件无关所以零维模型用c*加推力系数C_F的组合非常合适。推力输出也走同样的简化路径F C_F * A_t * pc。C_F由喷管面积比、燃气比热比和环境背压决定设计良好的喷管C_F通常在1.4到1.7左右。对方案阶段来说先按经验值取一个后面做喷管详细设计时再替换成随面积比变化的计算函数就行。3. 代码实现从参数表到ODE求解3.1 先准备好参数文件Matlab里我用结构体P装全部参数这样做的好处是后续做参数扫描时只需要循环修改结构体的字段函数内部不用动。参数写在一个setup脚本里如下% 推进剂与热力学参数 P.rho_p 1800; % 推进剂密度kg/m^3 P.a 5e-5; % 燃速系数m/s / Pa^n P.n 0.4; % 燃速压强指数 P.Tc 2800; % 绝热燃烧温度K P.R_gas 320; % 燃气气体常数J/(kg.K) P.c_star 1600; % 特征速度m/s P.CF 1.55; % 推力系数 % 药柱几何参数 P.d_i 0.04; % 药柱初始内径m P.D_o 0.09; % 药柱外径m P.L_grain 0.5; % 药柱长度m % 喷管与初始条件 P.A_t pi/4 * 0.03^2; % 喉部面积m^2 P.V0 2e-4; % 初始燃烧室自由容积m^3 P.p_init 101325; % 点火初始压强Pa注意这里的a5e-5是我按SI单位直接给的示例值对应的燃速在5 MPa下大约是24 mm/s属于高燃速复合推进剂的量级。如果你从文献查到的参数单位不是SI先做换算。3.2 ODE右端函数把方程翻译成代码整个模拟器的核心就是这个ODE右端函数。状态量我选了燃烧室压强pc和已烧蚀厚度e两个。e随时间的变化率就是燃速pc的变化率由上一节的微分方程给出function dydt srMotorODE(~, y, P) pc y(1); e y(2); d min(P.d_i 2*e, P.D_o); % 当前燃面直径 Ab pi * d * P.L_grain; % 当前燃面面积 rb P.a * pc^P.n; % Vieille燃速 Vc P.V0 pi/4 * (d^2 - P.d_i^2) * P.L_grain; dVdt pi/2 * d * P.L_grain * rb; % 自由容积变化率 m_dot_burn P.rho_p * Ab * rb; % 燃气生成率 m_dot_nozzle pc * P.A_t / P.c_star; % 喷管排出率 dpdt (P.R_gas * P.Tc / Vc) * (m_dot_burn - m_dot_nozzle) - (pc / Vc) * dVdt; dedt rb; dydt [dpdt; dedt]; end这里min函数是防止燃面直径超过药柱外径D_o避免ODE在燃尽后把几何尺寸推到物理上不可能的区间。3.3 事件检测燃尽如何停如果不做处理ODE求解器会在药柱燃尽后继续积分而此时d已经变成D_o燃面不再变生成为零数值解虽然不会崩但e还会被r_b持续推着走物理上就错了。正确做法是定义事件函数检测到已烧蚀厚度达到最大web厚度时终止积分function [value, isterminal, direction] burnEvent(~, y, P) e_max (P.D_o - P.d_i) / 2; value y(2) - e_max; isterminal 1; direction 0; end如果要做燃尽后的拖尾段也有办法在事件触发后把当前压强作为初值重新积分解一个纯排气方程dp/dt -(R_gas * Tc / V_final) * (p * A_t / c*)V_final取燃尽时刻的自由容积。我在实际项目中就是这样处理拖尾段的结果曲线和试车数据的尾部趋势对得上。3.4 后处理与绘图输出主求解段和后处理很直接opts odeset(RelTol, 1e-6, AbsTol, [1e4, 1e-8], Events, burnEvent); [t, y] ode15s((t,y) srMotorODE(t,y,P), [0 30], [P.p_init, 0], opts); pc y(:,1); e y(:,2); F P.CF * P.A_t * pc; % 推力曲线 figure; yyaxis left; plot(t, pc/1e6); ylabel(燃烧室压强 (MPa)); yyaxis right; plot(t, F/1000); ylabel(推力 (kN)); xlabel(时间 (s)); grid on;4. 数值求解的坑与验证从ode45到ode15s4.1 为什么我最后换了ode15s第一版代码我图省事直接用ode45跑结果点火上升段步长被压缩到微秒量级整个求解慢得离谱后面还出现过收敛警告。原因并不神秘点火瞬间燃烧室自由容积很小燃面很大燃气生成率瞬间起来而喷管排出项同时又强烈依赖压强两个时间尺度差了好几个数量级方程组呈现出明显的刚性特征。解决办法就是把求解器换成ode15s配置如下opts odeset(RelTol, 1e-6, AbsTol, [1e4, 1e-8], Events, burnEvent); [t, y] ode15s((t,y) srMotorODE(t,y,P), [0 30], [P.p_init, 0], opts);RelTol我设1e-6AbsTol里面压强分量设1e4 Pa厚度分量设1e-8 m。压强的绝对误差容限不能设太小因为量纲上压强的绝对量级是几兆帕1e4 Pa对应0.01 MPa的精度已经完全够用厚度e是微小量级所以给了更严格的1e-8。如果AbsTol设置不合适求解器会在某些很长的工作段反复取点拖慢速度。换完ode15s之后整个仿真在普通笔记本上秒级完成点火段的压强爬升过程也能清晰看到。这个经验后来也延续到了其它燃烧仿真项目里凡是涉及快速建立压强的动态过程我第一反应都是上隐式求解器。4.2 初始压强不能设成0点火初始条件有讲究运行中遇到的第一个报错竟然是压强不点火。我一开始把初始压强设成0结果Vieille公式里0的任意正指数次方都是0燃速为0燃烧生成项永远起不来仿真直接趴窝。这个坑写进代码注释里警示自己初始压强必须给一个能触发燃速起始的数值工程上一般直接给大气压101325 Pa代表点火药已经把燃烧室环境从真空或常压建立起来。更微妙的是初始自由容积V0。V0越小同样的燃气生成量对应的压强建立速度越快点火尖峰越高。如果V0取得过小仿真里会出现一个比稳态压强高出一倍多的尖锐峰值这在现实中对应点火冲击。反过来V0取得过大压强爬升变缓点火延迟感增强。做方案对比时我会把V0当作一个敏感参数单独扫描看它对待测发动机的峰值压强影响有多大。4.3 用解析解做自检别让曲线骗了你仿真跑通之后第一件事不是画图而是自检。零维内弹道有个很经典的解析解——稳态平衡压强p_bal ( (ρ_p * a * Ab * c*) / A_t )^(1/(1-n))这个公式在令dp/dt等于零、忽略dVdt项时可以得到。我通常用它验证仿真末段的稳定工作压强误差应在几个百分点以内。注意Ab是燃面面积由于内孔药柱是增面燃烧实际平衡压强会随着Ab增大而缓慢爬升所以严格说是一条缓慢上扬的平台而不是绝对水平线。自检时把初始Ab代入算一个p_bal再拿末端Ab代入算一个p_bal仿真曲线应当落在两者之间。另一条守恒关系是总冲。用trapz对F曲线做时间积分得到总冲再和推进剂质量乘以预计比冲对比I_total trapz(t, F); m_prop P.rho_p * pi/4 * (P.D_o^2 - P.d_i^2) * P.L_grain; Isp_avg I_total / (m_prop * 9.81);如果Isp_avg和推进剂理论比冲差超过5%就要回头查参数。我遇到过的情况是C_F给得过高导致推力虚高但压强曲线又正常这种问题不靠总冲守恒根本发现不了。5. 让模拟器更进一步点火瞬态、参数打靶与动态可视化5.1 点火瞬态模型基础版本里我把压强初始值设成大气压相当于是点火药已经把压强建立起来之后再交给主装药。要模拟完整的点火过程还得把点火药的质量生成率加进方程。常规做法是在ODE右端函数里增加一个点火药的燃烧生成项比如设定点火药在0到5毫秒内线性烧完产生一定质量的燃气等主装药压强达到着火阈值后主燃面才开始按Vieille公式产生燃气。这里有个很实用的小技巧用smoothstep或者分段线性函数来近似点火药的生成曲线比用阶跃更符合实际也更容易让ode15s稳定通过点火段。我经过对比发现点火药质量占推进剂总质量的0.1%到0.5%时点火峰值的相对量级和试车数据比较接近。5.2 蒙特卡洛参数打靶内弹道模型里几个参数天然具有散布性燃速系数a、压强指数n、喉部面积A_t。制造公差和工作环境的差异都会让这些参数偏离名义值。参数打靶的目的就是看偏差在合理范围内时最大压强和总冲的散布有多大是否还在结构裕度内。实现方式很简单写一个循环对a、n、A_t分别施加正态分布扰动比如标准差取名义值的2%到3%每次重新跑一遍仿真记录最大压强和总冲最后用直方图看分布。我在这个模拟器上加了这个功能之后很多方案评审问题可以直接用数据回答比如“喉部面积加工偏差3%的情况下最大压强有没有超过结构强度余量”。5.3 药柱烧蚀动态可视化最后一个让模拟器从工具变成演示程序的功能是药柱截面烧蚀过程的可视化。原理很简单每一帧都绘制当前的内孔圆和外圆填充推进剂区域颜色逐渐消失代表烧掉的部分。Matlab里用rectangle加faceColor配合循环帧更新十来行代码就能完成。我常用这个动画来做方案汇报和教学演示因为它能把抽象的内弹道曲线和药柱几何直观对应起来——当压强曲线出现明显爬升时画面里能看到内孔均匀扩大、燃面增大如果某一段压强表现异常也能从动画里快速判断是不是几何设计的问题。如果不想用动画也可以直接输出燃面面积随时间的变化曲线一样能说明问题。后续想和六自由度弹道模型耦合的话把推力曲线导出成CSV或.mat文件就行这个模拟器的输出格式我一开始就按这个接口预留了。本文还有配套的精品资源点击获取
简介面向固体火箭发动机设计与仿真领域的Matlab模拟器压缩包适合航天、机械及计算机仿真方向的研究者与工程师用于在物理试验前完成点火、燃烧、推进等环节的虚拟验证。包内共7个文件包含3个Matlab源程序覆盖主流程、发动机结构与推进剂模块2个br格式推进剂数据文件提供不同颗粒度燃烧参数另有说明文档与Git配置文件整体仅5KB结构精简。已有69人学习可用于快速搭建发动机几何模型、分析燃烧室与喷嘴流场、建立火药燃烧模型并求解运动方程涵盖从点火到推进的全过程仿真。借助Matlab可视化能力可直观观察温度、压力和流场变化支持参数优化与敏感性分析为设计验证和故障分析提供支撑对缩短研发周期、提升设计可靠性具有实用价值。 做固体火箭发动机内弹道仿真这件事我一开始并没有打算自己造轮子。当时手头有个小型固体火箭发动机的预研需求需要快速得到一组燃烧室压强曲线和推力曲线用来评估装药方案到底靠不靠谱。翻了一圈现成工具后我发现一个尴尬局面通用内弹道程序要么药型库固定想改几何参数得动源代码要么干脆是个黑盒算完给你一张曲线图中间压强怎么建立、燃面怎么发展、自由容积怎么变化全都看不见。用CFD又太重网格、湍流模型、两相流一轮下来项目周期根本扛不住。最后我决定在Matlab里手写一个固体火箭发动机模拟器基于零维内弹道模型把燃速模型、药柱几何、喷管流率和ODE求解完整串起来。跑通后十来分钟就能出一组结果而且每一步物理过程都能拆开看后来这个模拟器陪我改了好几轮方案也提前筛掉过两个参数上明显不合理的药型。这篇内容适合对固体火箭发动机有基本概念、想真正把零维内弹道代码跑通的工程师或学生我会把物理模型、代码实现、数值坑和验证方法一次讲清楚。1. 为什么自己写一个Matlab模拟器不直接套现成软件1.1 这个模拟器能算什么固体火箭发动机内弹道仿真的核心产出其实就是两条曲线燃烧室压强随时间的变化以及推力随时间的变化。这两条曲线几乎决定了发动机的所有主要性能——总冲、平均推力、最大压强、工作时间、装药是否合理。再往下还能导出燃面面积随时间的退化、推进剂剩余质量、点火瞬间的压强冲击高度。对预研阶段来说这个量级的信息已经完全够用。你不需要知道燃烧室内部哪里的流场有回流、哪里的温度更高因为零维模型本身就是把整个燃烧室当作一个充分混合的控制体压强在任意瞬间处处相等。这个假设听起来很粗暴但对绝大多数固体火箭发动机的设计迭代来说它给出的结果已经相当能打。用这套模拟器你可以做的事情包括比较不同药柱内径下的压强平台段评估喉部面积变大后对工作压强的影响估算点火瞬态的峰值压强计算总冲和平均比冲。这些都是在方案阶段必须回答的问题而且这些问题用手算只有稳态解看不到动态过程用CFD又像是在用大炮打蚊子。1.2 为什么选Matlab而不是专用软件有人会问NASA的CEA、商业软件里的固体发动机模块都能算为什么要自己在Matlab里折腾。我当时的判断很简单对比项Matlab模拟器专用内弹道程序通用CFD学习成本低核心方程能自己推导中需要熟悉输入格式高需要网格和模型经验可修改性完全开放改一行代码就换一种药型受程序架构限制改几何就要重新建模计算速度秒级到分钟级快小时级到天级结果透明度每个中间变量都能看只有最终曲线后处理复杂前期投入只花时间可能涉及授权人力物力都大Matlab还有个天然优势ODE求解器非常成熟数值刚性检测、事件触发、误差控制都有完整方案而且绘图、数据处理、参数扫描都在同一个环境里完成不用在Python、Origin、Excel之间来回倒腾。如果你在高校或研究所Matlab基本是标配不存在环境门槛。我把这套模拟器写成纯脚本加函数不依赖Simulink装个基础版Matlab就能跑。2. 零维内弹道模型的物理方程从质量守恒到压强微分2.1 核心闭环三个变量决定了燃烧室压强零维内弹道模型的根基就是燃烧室内的质量守恒。压强随时间的变化本质上由三个量在博弈推进剂燃烧产生的气体、喷管排出的气体、以及燃烧室自由容积的变化。写成微分方程就是dp_c/dt (R_gas * Tc / Vc) * (ṁ_gen - ṁ_nozzle) - (pc / Vc) * dVc/dt这里每一项都要拆开看。第一项里的R_gas是燃气的气体常数Tc是燃烧温度Vc是当前燃烧室自由容积。ṁ_gen是燃烧生成的气体质量流率正比于推进剂密度、当前燃面面积和燃速ṁ_nozzle是喷管排出的质量流率。最后一项里的dVc/dt是自由容积的变化率因为药柱不断烧掉内孔扩大留给气体的空间在不断变大。这个方程要说直观也直观生成比排出多压强升高自由容积增大压强倾向下降。两者达到平衡时dp_c/dt为零发动机就进入稳态工作段。这里我强调一点很多简化模型会忽略dVc/dt这一项但点火段和薄药柱的末段这一项的影响会被放大如果完全不考虑压强曲线会显得过于平直和实测对不上。2.2 燃速模型与药柱几何参数全在几何里燃速采用固体推进剂最经典的Vieille经验公式r_b a * pc^na是燃速系数n是压强指数。压强指数n的大小直接影响发动机的稳定性一般复合推进剂在0.3到0.5之间。指数越高压强波动越容易被放大所以设计时通常想办法压低n。参数单位是最大的坑。很多资料里给出a的单位是mm/s·MPa^-n但SI体系下压强单位是Pa直接代进去结果差好几个数量级。换算方法很简单如果查到的a_cmps是以mm/s和MPa为基准的转成SI要乘上(1e6)^(-n)因为1 MPa等于1e6 Pa。比如a1.2 mm/s·MPa^-0.4转成SI就是1.2e-3再乘以1e6的负0.4次方。这个换算我见过太多人栽跟头算出来的压强要么离谱地高要么离谱地低。药柱几何我以最常用的圆柱内孔药柱为例两端做阻燃包覆只有内孔表面参与燃烧。当前燃面直径d等于初始内径d_i加上两倍的已烧蚀厚度e燃面面积Ab π * d * L_grain。随着燃烧推进d变大Ab也跟着变大所以燃面是增面燃烧平衡压强会缓慢上升直到外层烧完。自由容积同步更新Vc V0 π/4 * (d^2 - d_i^2) * L_grainV0是装药前就存在的初始空腔容积也包括点火器空间。这个V0对点火瞬态影响极大后面我会专门说。2.3 喷管流率与特征速度c*喷管排出项写的是ṁ_nozzle pc * A_t / c*其中A_t是喉部面积c*是特征速度。这里没有真的去解喷管内流动而是用c*这个参数把燃烧室到喉部之间的能量转化打包描述。c*可以由燃烧温度、燃气常数、比热比算出c* sqrt(R_gas * Tc / γ) / sqrt((2/(γ1))^((γ1)/(γ-1)))工程上更常见的是直接用推进剂手册里的实验值比如典型的AP/HTPB复合推进剂c*大约在1500到1650 m/s。c*只反映燃烧气体的能量水平跟喷管扩张比和出口条件无关所以零维模型用c*加推力系数C_F的组合非常合适。推力输出也走同样的简化路径F C_F * A_t * pc。C_F由喷管面积比、燃气比热比和环境背压决定设计良好的喷管C_F通常在1.4到1.7左右。对方案阶段来说先按经验值取一个后面做喷管详细设计时再替换成随面积比变化的计算函数就行。3. 代码实现从参数表到ODE求解3.1 先准备好参数文件Matlab里我用结构体P装全部参数这样做的好处是后续做参数扫描时只需要循环修改结构体的字段函数内部不用动。参数写在一个setup脚本里如下% 推进剂与热力学参数 P.rho_p 1800; % 推进剂密度kg/m^3 P.a 5e-5; % 燃速系数m/s / Pa^n P.n 0.4; % 燃速压强指数 P.Tc 2800; % 绝热燃烧温度K P.R_gas 320; % 燃气气体常数J/(kg.K) P.c_star 1600; % 特征速度m/s P.CF 1.55; % 推力系数 % 药柱几何参数 P.d_i 0.04; % 药柱初始内径m P.D_o 0.09; % 药柱外径m P.L_grain 0.5; % 药柱长度m % 喷管与初始条件 P.A_t pi/4 * 0.03^2; % 喉部面积m^2 P.V0 2e-4; % 初始燃烧室自由容积m^3 P.p_init 101325; % 点火初始压强Pa注意这里的a5e-5是我按SI单位直接给的示例值对应的燃速在5 MPa下大约是24 mm/s属于高燃速复合推进剂的量级。如果你从文献查到的参数单位不是SI先做换算。3.2 ODE右端函数把方程翻译成代码整个模拟器的核心就是这个ODE右端函数。状态量我选了燃烧室压强pc和已烧蚀厚度e两个。e随时间的变化率就是燃速pc的变化率由上一节的微分方程给出function dydt srMotorODE(~, y, P) pc y(1); e y(2); d min(P.d_i 2*e, P.D_o); % 当前燃面直径 Ab pi * d * P.L_grain; % 当前燃面面积 rb P.a * pc^P.n; % Vieille燃速 Vc P.V0 pi/4 * (d^2 - P.d_i^2) * P.L_grain; dVdt pi/2 * d * P.L_grain * rb; % 自由容积变化率 m_dot_burn P.rho_p * Ab * rb; % 燃气生成率 m_dot_nozzle pc * P.A_t / P.c_star; % 喷管排出率 dpdt (P.R_gas * P.Tc / Vc) * (m_dot_burn - m_dot_nozzle) - (pc / Vc) * dVdt; dedt rb; dydt [dpdt; dedt]; end这里min函数是防止燃面直径超过药柱外径D_o避免ODE在燃尽后把几何尺寸推到物理上不可能的区间。3.3 事件检测燃尽如何停如果不做处理ODE求解器会在药柱燃尽后继续积分而此时d已经变成D_o燃面不再变生成为零数值解虽然不会崩但e还会被r_b持续推着走物理上就错了。正确做法是定义事件函数检测到已烧蚀厚度达到最大web厚度时终止积分function [value, isterminal, direction] burnEvent(~, y, P) e_max (P.D_o - P.d_i) / 2; value y(2) - e_max; isterminal 1; direction 0; end如果要做燃尽后的拖尾段也有办法在事件触发后把当前压强作为初值重新积分解一个纯排气方程dp/dt -(R_gas * Tc / V_final) * (p * A_t / c*)V_final取燃尽时刻的自由容积。我在实际项目中就是这样处理拖尾段的结果曲线和试车数据的尾部趋势对得上。3.4 后处理与绘图输出主求解段和后处理很直接opts odeset(RelTol, 1e-6, AbsTol, [1e4, 1e-8], Events, burnEvent); [t, y] ode15s((t,y) srMotorODE(t,y,P), [0 30], [P.p_init, 0], opts); pc y(:,1); e y(:,2); F P.CF * P.A_t * pc; % 推力曲线 figure; yyaxis left; plot(t, pc/1e6); ylabel(燃烧室压强 (MPa)); yyaxis right; plot(t, F/1000); ylabel(推力 (kN)); xlabel(时间 (s)); grid on;4. 数值求解的坑与验证从ode45到ode15s4.1 为什么我最后换了ode15s第一版代码我图省事直接用ode45跑结果点火上升段步长被压缩到微秒量级整个求解慢得离谱后面还出现过收敛警告。原因并不神秘点火瞬间燃烧室自由容积很小燃面很大燃气生成率瞬间起来而喷管排出项同时又强烈依赖压强两个时间尺度差了好几个数量级方程组呈现出明显的刚性特征。解决办法就是把求解器换成ode15s配置如下opts odeset(RelTol, 1e-6, AbsTol, [1e4, 1e-8], Events, burnEvent); [t, y] ode15s((t,y) srMotorODE(t,y,P), [0 30], [P.p_init, 0], opts);RelTol我设1e-6AbsTol里面压强分量设1e4 Pa厚度分量设1e-8 m。压强的绝对误差容限不能设太小因为量纲上压强的绝对量级是几兆帕1e4 Pa对应0.01 MPa的精度已经完全够用厚度e是微小量级所以给了更严格的1e-8。如果AbsTol设置不合适求解器会在某些很长的工作段反复取点拖慢速度。换完ode15s之后整个仿真在普通笔记本上秒级完成点火段的压强爬升过程也能清晰看到。这个经验后来也延续到了其它燃烧仿真项目里凡是涉及快速建立压强的动态过程我第一反应都是上隐式求解器。4.2 初始压强不能设成0点火初始条件有讲究运行中遇到的第一个报错竟然是压强不点火。我一开始把初始压强设成0结果Vieille公式里0的任意正指数次方都是0燃速为0燃烧生成项永远起不来仿真直接趴窝。这个坑写进代码注释里警示自己初始压强必须给一个能触发燃速起始的数值工程上一般直接给大气压101325 Pa代表点火药已经把燃烧室环境从真空或常压建立起来。更微妙的是初始自由容积V0。V0越小同样的燃气生成量对应的压强建立速度越快点火尖峰越高。如果V0取得过小仿真里会出现一个比稳态压强高出一倍多的尖锐峰值这在现实中对应点火冲击。反过来V0取得过大压强爬升变缓点火延迟感增强。做方案对比时我会把V0当作一个敏感参数单独扫描看它对待测发动机的峰值压强影响有多大。4.3 用解析解做自检别让曲线骗了你仿真跑通之后第一件事不是画图而是自检。零维内弹道有个很经典的解析解——稳态平衡压强p_bal ( (ρ_p * a * Ab * c*) / A_t )^(1/(1-n))这个公式在令dp/dt等于零、忽略dVdt项时可以得到。我通常用它验证仿真末段的稳定工作压强误差应在几个百分点以内。注意Ab是燃面面积由于内孔药柱是增面燃烧实际平衡压强会随着Ab增大而缓慢爬升所以严格说是一条缓慢上扬的平台而不是绝对水平线。自检时把初始Ab代入算一个p_bal再拿末端Ab代入算一个p_bal仿真曲线应当落在两者之间。另一条守恒关系是总冲。用trapz对F曲线做时间积分得到总冲再和推进剂质量乘以预计比冲对比I_total trapz(t, F); m_prop P.rho_p * pi/4 * (P.D_o^2 - P.d_i^2) * P.L_grain; Isp_avg I_total / (m_prop * 9.81);如果Isp_avg和推进剂理论比冲差超过5%就要回头查参数。我遇到过的情况是C_F给得过高导致推力虚高但压强曲线又正常这种问题不靠总冲守恒根本发现不了。5. 让模拟器更进一步点火瞬态、参数打靶与动态可视化5.1 点火瞬态模型基础版本里我把压强初始值设成大气压相当于是点火药已经把压强建立起来之后再交给主装药。要模拟完整的点火过程还得把点火药的质量生成率加进方程。常规做法是在ODE右端函数里增加一个点火药的燃烧生成项比如设定点火药在0到5毫秒内线性烧完产生一定质量的燃气等主装药压强达到着火阈值后主燃面才开始按Vieille公式产生燃气。这里有个很实用的小技巧用smoothstep或者分段线性函数来近似点火药的生成曲线比用阶跃更符合实际也更容易让ode15s稳定通过点火段。我经过对比发现点火药质量占推进剂总质量的0.1%到0.5%时点火峰值的相对量级和试车数据比较接近。5.2 蒙特卡洛参数打靶内弹道模型里几个参数天然具有散布性燃速系数a、压强指数n、喉部面积A_t。制造公差和工作环境的差异都会让这些参数偏离名义值。参数打靶的目的就是看偏差在合理范围内时最大压强和总冲的散布有多大是否还在结构裕度内。实现方式很简单写一个循环对a、n、A_t分别施加正态分布扰动比如标准差取名义值的2%到3%每次重新跑一遍仿真记录最大压强和总冲最后用直方图看分布。我在这个模拟器上加了这个功能之后很多方案评审问题可以直接用数据回答比如“喉部面积加工偏差3%的情况下最大压强有没有超过结构强度余量”。5.3 药柱烧蚀动态可视化最后一个让模拟器从工具变成演示程序的功能是药柱截面烧蚀过程的可视化。原理很简单每一帧都绘制当前的内孔圆和外圆填充推进剂区域颜色逐渐消失代表烧掉的部分。Matlab里用rectangle加faceColor配合循环帧更新十来行代码就能完成。我常用这个动画来做方案汇报和教学演示因为它能把抽象的内弹道曲线和药柱几何直观对应起来——当压强曲线出现明显爬升时画面里能看到内孔均匀扩大、燃面增大如果某一段压强表现异常也能从动画里快速判断是不是几何设计的问题。如果不想用动画也可以直接输出燃面面积随时间的变化曲线一样能说明问题。后续想和六自由度弹道模型耦合的话把推力曲线导出成CSV或.mat文件就行这个模拟器的输出格式我一开始就按这个接口预留了。本文还有配套的精品资源点击获取
