2-57 气缸从充气到放气，压力、温度和流量会怎样随时间变化？

MATLAB ODE45 / ODE23t 气缸充放气热力学与无因次仿真程序

一句话介绍： 这套程序以 30 L 定容气缸为对象，用分段质量流量模型描述气源与容腔之间的气体交换，再分别建立多变/等温/绝热充气、带传热充气以及无因次充气和放气微分方程，最终得到压力、温度和归一化流量随过程推进的动态变化。

先看它解决什么问题

气缸充气和放气看起来只是“气体进来、气体出去”，但真正做动态仿真时，至少有三件事会同时变化：

容腔压力不断改变，入口与出口压比随之变化；

质量流量不是常数，会经历临界流动与非临界流动区间；

气体压缩或膨胀会改变温度，而壁面传热又会把温度拉回环境温度。

当前工程把这些因素拆成四条可独立运行的仿真链：

供气参数与初始状态 → 分段质量流量 → 压力/温度状态方程 → ODE 数值积分 → 压力、温度、流量动态曲线

其中 solve1SYS.m 用于多变过程及等温/绝热对比，solve2SYS.m 加入传热项，solve3SYS.m 和 solve4SYS.m 则分别处理无因次充气与放气。

1. 当前工程里真正有哪些仿真文件？

压缩包没有统一的总入口，而是由四个独立 MATLAB 主函数组成，每个 .m 文件同时包含主函数和对应的局部微分方程函数。

solve4SYS.asv 是 MATLAB 自动保存的备份文件，不参与正常运行。

这意味着它不是一个 GUI，也不依赖外部数据文件；四个主函数就是四个直接的实验入口。

2. 当前程序模拟的气缸和工况是什么？

源码中最核心的对象是一个 30 L 定容气缸。主要参数为：

气体常数 R = 287 J/(kg·K)；

容积 V = 30 L；

声速流导 C = 0.2 dm³/(s·bar)；

临界压力比 b = 0.4；

初始/环境压力 P₀ = 1 bar；

当前源码供气压力 p₁ 或 Ps = 6 bar；

供气温度、环境温度均取 293 K；

标准状态空气密度采用 1.185×10⁻³ kg/dm³。

源码注释还给出有效截面积与声速流导的换算关系：Se ≈ 5.03C，因此当前 C=0.2 对应约 1.006 mm²。

这里有一个需要按资料优先级处理的差异：配套 Word 题设中出现了 Ps = 0.5 MPa，而当前四个 MATLAB 源文件实际使用的是 6 bar。本文以压缩包中的当前源码为准，因此后续按 6 bar 讨论。

3. 气体进入气缸时，质量流量为什么不是常数？

当前程序先根据容腔压力与气源压力的比值决定质量流量。令压力比为：

r = p/p₁

最大质量流量按源码可写成：

Gₘₐₓ = ρ₀ C p₁ √(293/θ₁)

当前 θ₁ = 293 K，所以根号项为 1。

程序随后按临界压力比 b=0.4 分段：

当 r ≤ b 时，G = Gₘₐₓ；

当 b < r ≤ 1 时，G = Gₘₐₓ × √[1 - ((r - b)/(1 - b))²]；

当容腔压力达到或超过气源压力后，理想情况下质量流量应降为 0。

直观上就是：气缸刚开始压力低，供气能力基本保持最大；随着气缸压力升高，压差变小，流量逐渐下降；两端压力接近后，充气过程结束。

这条分段流量关系是后面四个模型共同的基础。

4. 等温、绝热和多变充气在程序里是怎样统一起来的？

solve1SYS.m 用一个指数 k 把不同热力过程放进同一组状态方程。当前默认：

k = 1.4 —— 源码标记为绝热过程。

如果把这一行改成：

k = 1 —— 就得到等温过程。

当前微分方程为：

dp/dt = kGRθ/V

dθ/dt = (k - 1)GRθ²/(pV)

这里最关键的不是公式本身，而是 k 对温度方程的影响。

当 k=1 时：

dθ/dt = 0

所以温度始终保持初值 293 K，压力只随进入气缸的质量增加而升高。

当 k=1.4 时，第二个方程不再为零，气体在充气过程中升温。把两个状态方程相除，还能得到当前代码对应的压力—温度关系：

θ/θ₀ = (p/p₀)^[(k - 1)/k]

因此这套程序所谓的“多变过程”，本质上就是通过不同 k 值改变充气过程中的热状态响应。

5. 默认绝热与等温充气会得到什么差别？

按 solve1SYS.m 当前方程和默认参数做等价数值复算，初始状态均为：

p₀ = 1 bar，θ₀ = 293 K

当默认使用 k=1.4 时，压力最终趋近当前供气压力 6 bar，温度随压力上升，最终约为：

θ ≈ 488.87 K

如果把 k 改为 1，压力同样趋近 6 bar，但温度保持：

θ = 293 K

按当前方程复算，压力接近 5.999 bar 的时间大约为：

这个差别很直观：在当前模型中，绝热充气时温度升高会同时进入压力变化方程，因此压力上升得更快；等温模型把所有温度变化直接压制为零。

6. 如果气缸会和环境换热，模型怎样变化？

solve2SYS.m 在质量与能量平衡中加入了传热项。

程序先计算定容、定压比热：

Cv = R/(K - 1)

Cp = Cv + R

再用等效表面积 Sh 和环境温差构造壁面换热：

Q̇ = hSh(θa - θ)

于是压力变化不再只由进气决定，而是同时受到进气焓流和壁面传热影响；温度变化也同时包含：

进入气体带来的能量 → 容腔内气体内能变化 → 与环境之间的热交换

这正是 solve2SYS.m 与前一个模型最重要的区别：它不再把过程限定成理想等温或理想绝热，而是让温度根据能量平衡自己变化。

配套题设写的是换热系数 20 W/(m²·K)，当前 solve2SYS.m 中变量写成 h=20e3，而无因次程序中又使用 h=20。因此这几个文件的数值表示并没有完全统一，迁移到新的物理工况时不应只机械照搬 h，而应先统一面积、压力、体积和热量的单位体系。

7. 为什么程序还要再做一套无因次充气模型？

solve3SYS.m 不再直接用 bar 和 K 作为主要状态，而是把压力、温度和流量归一化。

当前初始压力比为：

P₀* = P₀/Ps = 1/6

初始无因次温度为：

θ₀* = 1

充气时，无因次流量 G* 仍采用 b=0.4 的分段形式：低压阶段保持 1，随后随压力比增加而逐渐下降到 0。

当前状态方程为：

dP*/dt = K G* + Ka(1 - θ*)

dθ*/dt = (θ*/P*)[(K - θ*)G* + Ka(1 - θ*)]

其中 Ka 用来表示传热时间尺度相对于充气时间尺度的影响。代码虽然计算了 Tp 和 Th，并在注释中写有 Ka = Tp/Th ≈ 1.0846，但真正参与当前仿真的值是：

Ka = 10

所以当前无因次结果反映的是 Ka=10 的工况，而不是注释中的 1.0846。

8. 放气模型和充气模型最大的不同在哪里？

solve4SYS.m 从高压状态开始：

P*(0) = 1，θ*(0) = 1

环境压力仍对应：

P₀* = 1/6

放气时真正决定出口流量的是外界压力与容腔压力之比：

r = P₀*/P*

因此 G* 不仅包含临界压力比修正，还乘有当前压力和温度项。当前程序的核心形式可以概括为：

G* = P*/√θ* × 流量修正项

随后压力方程中的流出项带负号：

dP*/dt = -K G*θ* + Ka(1 - θ*)

温度方程也对应气体流出造成的膨胀降温，再叠加环境传热恢复项。

所以充气与放气的物理方向正好相反：

充气：压力升高、温度先上升、流量逐渐衰减

放气：压力降低、温度先下降、流量逐渐衰减

9. 当前无因次模型能看到什么典型动态？

按 solve3SYS.m、solve4SYS.m 当前方程以及实际生效的 Ka=10 做等价复算，可以得到很清晰的趋势。

无因次充气

P* 从 1/6 上升到约 1；

G* 从最大流量逐渐降到 0；

θ* 先升高，峰值约 1.036，随后在传热作用下回到约 1。

若按环境温度 293 K 换算，这个峰值约对应 303.65 K。

无因次放气

P* 从 1 下降到约 1/6；

G* 最终同样衰减到 0；

θ* 先下降，最低约 0.970，随后回到约 1。

按 293 K 换算，最低温度约为 284.14 K。

这组结果说明 Ka=10 时壁面换热比较强，温度偏离环境值后会较快恢复，而压力仍完成从初始状态到目标压力状态的过渡。

需要注意：solve3SYS.m 和 solve4SYS.m 明确写的是“无因次仿真”，代码虽然计算了时间尺度 Tp，但没有把积分变量显式乘回物理时间；横坐标仍标成 t/s。因此这两组曲线更适合用于比较无因次动态趋势，不宜在未补充时间尺度换算前直接把横坐标解释成严格的实际秒数。

10. 这套程序适合用在哪里？

气动系统充放气过程教学

程序把质量流量分段、压力变化、温度变化和换热直接放进微分方程，非常适合用来理解“为什么充气流量会下降”“为什么绝热过程会升温”“为什么放气会降温”等基础动态现象。

等温、绝热与传热模型对比

solve1SYS.m 可以直接通过 k=1 和 k=1.4 对比两个极限假设，solve2SYS.m 又进一步加入热交换，适合做不同热力假设下的响应比较。

气缸和气罐参数敏感性分析

容积、供气压力、流导、临界压力比和传热强度都集中写在源码开头，适合做离线参数扫描，观察充气速度、放气速度和温度峰值如何改变。

气动瞬态模型的前期验证

在建立更复杂的阀—管路—气缸系统之前，可以先用这套定容模型检查基本的质量守恒、能量变化和压比流量规律。不过当前工程属于 MATLAB 数值仿真脚本，不是实时控制程序，也没有阀芯动力学、管路压力损失或执行机构机械负载模型。

11. 哪些参数最值得修改？

当前程序里真正值得主动做实验的，是下面这些会直接改变压力、温度或流量动态的参数。

其中 Ka 最需要特别注意：源码已经计算 Tp/Th，但当前没有使用该计算结果，而是直接固定为 10。做理论复现或参数对比时，应明确自己到底要使用当前程序值，还是使用 Tp/Th 得到的值。

12. 怎么运行这套程序？

把压缩包中的 MATLAB 文件放在同一工作目录，并将该目录设为 MATLAB 当前路径后，可以分别运行：

solve1SYS、solve2SYS、solve3SYS、solve4SYS

其中：

运行 solve1SYS 时，默认是 k=1.4；要看等温结果，把 k=1.4 注释掉并启用 k=1；

solve3SYS 和 solve4SYS 当前真正使用的是 Ka=10；如果要复现注释中的其他传热参数，需要手动修改该行；

四个文件都已经把对应的微分方程函数写在同一个 .m 文件内部，不需要再额外寻找独立的 solve1.m、solve2.m、solve3.m 或 solve4.m；

程序使用 MATLAB 自带的 ode45 和 ode23t 数值求解器，没有外部数据文件依赖。

还有一个会影响 solve2SYS.m 稳定使用的源码细节：当前压缩包里的分段流量判断写成了 else 后再单独计算 (y(1)/p1<=1)，而压力超过气源压力时的 G=0 分支被注释掉了。按物理分段关系使用时，应把中间判断整理为真正的 elseif y(1)/p1 <= 1，并保留最后的 else → G=0，否则数值积分一旦越过供气压力，根号项可能进入无效区间。

13. 一句话看懂这个项目

这是一个定容气缸充放气动态仿真项目：它以 30 L 容腔、C=0.2 dm³/(s·bar)、b=0.4、初始压力 1 bar 和当前源码供气压力 6 bar 为基础，先通过压比建立分段质量流量，再用 ODE45/ODE23t 求解压力与温度状态；solve1SYS.m用k=1/1.4对比等温与绝热充气，solve2SYS.m加入壁面传热，solve3SYS.m与solve4SYS.m用Ka=10分别描述无因次充气和放气。按当前方程等价复算，默认绝热充气由 1 bar 升至 6 bar 时温度约升至 488.87 K，无因次充气温度峰值约为 1.036、放气最低约为 0.970；这套程序更适合用于气动热力过程教学、参数敏感性分析和气缸/气罐瞬态模型的前期验证。