稳态选型做完之后很多人以为自己已经把系统摸透了直到要做启停、变负荷、除霜、电子膨胀阀整定这几件事才发现手里那套基于设计工况的计算表格完全不够用。PyVCRSDyn 这个项目要解决的正是这一段空白用 Python 搭一套面向蒸汽压缩制冷与热泵循环的制冷系统瞬态仿真框架把压缩机、冷凝器、节流装置、蒸发器连成一个可积分的整体让压力、过热度、壁温这些量随时间跑起来。它不是又一个稳态选型工具也不是把现成商业软件的界面翻译一遍而是把瞬态建模里最麻烦的几件事——移动边界的段长切换、两相区物性、刚性方程求解、事件触发——用一套结构清晰、可以自己动手改的代码组织起来。适合谁来读手上正在做热泵控制策略验证、冷库温度波动分析、车用空调启停冲击评估或者单纯想搞明白为什么我的仿真跑三步就发散的人。1. 稳态选型算不出的那几秒才是这套框架存在的理由1.1 设计工况之外系统绝大部分时间都不在稳态做过换热器选型的人都熟悉那套流程确定制冷量、蒸发温度、冷凝温度查压焓图算质量流量配换热面积选压缩机排量。这套方法在选型阶段非常好用也几乎不会出错。问题在于它默认了一件事——系统有无限长的时间去达到那个平衡点。真实设备不是这样运行的。一台变频热泵在冬季供暖时压缩机转速可能每分钟都在变一台冷库在货物入库的半小时内蒸发器负荷会翻倍一台车用空调在怠速和高速之间的冷凝风量差异足以让高压侧压力在一个呼吸周期内上下摆动。这些过程里稳态计算给的只是一个终点而你真正关心的往往是从 A 点到 B 点中间发生了什么。压缩机刚启动的那几秒排气温度会不会冲到保护阈值电子膨胀阀从 30% 开到 60%过热度要多久才回落中间会不会先掉到负值导致回液除霜结束时换热器里积的那点液态制冷剂被压缩机吸走会造成多大的液击风险。这些问题没有一个能靠稳态表格回答。我见过不少团队的做法是先用稳态软件算出设计点然后拿一个简化的集总模型凑合一下把整个换热器当成一个温度均匀的大热容。这种模型在负荷缓慢变化时勉强能用但只要涉及启停或者阀门动作误差就会大到没法看——因为它完全丢掉了换热器内部的沿程分布而过热度这种东西恰恰是由分布决定的。1.2 三个时间尺度叠在一起决定了你必须用隐式求解瞬态仿真能不能跑起来很大程度上取决于你有没有认真对待系统里的时间尺度分布。制冷循环里至少有三个量级完全不同的惯性环节主导物理量典型时间常数压缩机转速与阀开度响应机械/电气惯性、控制器周期毫秒到 1 秒换热器管壁与制冷剂存量金属热容、两相区质量1 秒到 30 秒房间、箱体、货物、蓄冷介质大热容对象数分钟到数小时粗算一下管壁这一项一段 0.5 kg 的铜管比热约 385 J/(kg·K)热容就是 190 J/K 左右如果这一段的等效换热系数乘面积在 200 W/K 量级时间常数就是 1 秒上下。而一个几立方米的冷库箱体加货物的等效热容轻松到 10^5 J/K 量级时间常数上千秒。两端的比值达到 10^3 到 10^4这就是典型的刚性系统。刚性系统的直接后果是显式积分方法的步长被最快的那个模态死死锁住。你想模拟一小时的除霜周期结果步长被限制在毫秒级要算几百万步才能走完。这不是算力问题是方法选错了。所以这套框架从第一行代码起就得按 DAE/刚性 ODE 来设计用隐式方法、给雅可比、处理事件而不是先写个欧拉法跑通了再说。我见过太多项目在后期被迫推倒重来原因就是早期选了显式积分等到想加管壁热容的时候发现完全跑不动。1.3 我把第一版的需求砍到只剩四条开源框架最容易犯的毛病是贪多。第一版什么都想支持最后什么都不能用。我给 PyVCRSDyn 的第一版定了一个非常克制的目标只有四条能跑单级蒸汽压缩循环一台压缩机、一个冷凝器、一个节流装置、一个蒸发器闭式连接工质单组分。换热器支持移动边界法蒸发器两段两相、过热冷凝器三段过热、两相、过冷段长作为状态量随时间变化。状态量固定为压力和比焓这是压焓图思维的直接延续也是物性反查最稳的组合。有事件机制和一致性初始化能从稳态解起步能在段长趋零时切换模式而不崩。这四条之外的比如多回路、混合工质、润滑油、结霜、风扇动态全部留接口不做实现。理由很简单一个能被验证的最小闭环价值远大于一个功能齐全但没人敢信的黑箱。你把这四条做扎实了后面加东西是增量工作做不扎实后面每加一个功能都是在补窟窿。2. 部件、端口与状态量先把骨架立对2.1 每个部件都是一个带守恒方程的端口盒子框架的骨架我参考了面向对象的系统建模思路每个物理部件是一个类对外只暴露两类东西——端口和状态量。端口负责传质量流量、比焓、压力这些流过边界的量状态量是部件内部真正被积分的东西。部件之间靠端口连接连接规则只有一条流入某端口的量等于流出对端口的量压力相等。这么设计的好处是可替换性。你想把移动边界法的蒸发器换成有限体积法的版本只要端口签名一致上层回路一行都不用改。想加一个回热器也只是在回路里插一个两头端口的新类。class Port: 端口只承载流过的量不承载状态 def __init__(self, name): self.name name self.mdot 0.0 # 质量流量 kg/s流入为正 self.h 0.0 # 比焓 J/kg self.p 0.0 # 压力 Pa class Component: 所有部件的基类 def __init__(self, name): self.name name self.inlet Port(in) self.outlet Port(out) def state_size(self): 返回本部件的状态量个数 raise NotImplementedError def initial_state(self, y0, offset): 把初始值写进全局状态向量 raise NotImplementedError def residuals(self, t, y, dy, ctx): 返回本部件贡献的方程残差DAE 形式 F(t,y,dy)0 raise NotImplementedError这里有个设计取舍值得说清楚我选择让部件返回残差而不是右端项。原因是 DAE 形式更通用——压缩机和节流阀本质上是代数方程没有自己的微分状态如果强行写成 dy/dt f(...) 的形式就得给它编一个假的容量项物理上不干净。用残差形式代数方程就是残差恒等于零微分方程就是残差等于 dy/dt 减去右端项两种统一处理。2.2 状态量怎么挑压力、焓、壁温、转速状态量的选择直接决定了这套框架好不好用、稳不稳。选压力加比焓这一组是我试过几种组合之后的结论。先说为什么不用温度加压力。两相区里温度和压力是绑死的你给了压力温度就定了反过来也一样。如果两个都当状态两相区里就出现了冗余雅可比矩阵直接奇异。这是新手最容易踩的坑报错信息通常是矩阵奇异或者步长被压到极小很难一眼看出是状态选重了。也不用密度加内能。数学上没问题但工程上不方便你从传感器拿到的、从选型表查到的、跟客户沟通的全是压力和温度、过热度。用焓做状态压焓图上一条竖线就是等压过程物理直觉非常顺。具体到各部件蒸发器状态是出口压力、两相段出口比焓、过热段出口比焓、两段壁温、两段长度。过热度不用单独当状态它就是出口焓减去饱和气焓再换算出来的东西。冷凝器状态是压力、三个段的出口比焓、三段壁温、三段长度。压缩机稳态版没有状态只有代数方程如果要模拟变频响应加一个转速的一阶惯性状态就一个。节流装置纯代数状态为零。开度作为外部输入或者控制器的输出。壁温必须当状态这一点不能省。很多简化模型直接把管壁热容折算进制冷剂里看起来省了几个状态实际上把制冷剂响应快、管壁响应慢这个关键的动态分离给抹掉了算出来的过热度响应会明显偏快。2.3 代数环放在哪一层解决定了你能不能跑起来闭式循环里天然存在代数环。蒸发压力影响压缩机吸气密度吸气密度影响质量流量质量流量影响蒸发器里的存量存量又反过来影响压力。如果你老老实实按压力是代数变量、由总质量守恒确定来写得到的就是一个纯 DAE指标可能大于 1求解器会非常难受。我在框架里用了两个手段来缓解。第一个是虚拟容积法给每个压力节点配一个很小的容量项把质量守恒的代数约束改写成压力的快速微分方程。$$ C_p \frac{dp}{dt} \sum \dot m_{in} - \sum \dot m_{out}, \qquad C_p V \cdot \left(\frac{\partial \rho}{\partial p}\right)_h $$这个容量系数本质上就是这点容积里压力变一点能存多少质量。它把代数环变成了一个时间常数极小的微分方程求解器只要把它当刚性模态处理就行比处理纯代数约束稳得多。代价是引入了一点非物理的快速动态所以容量系数不能取得太小一般让对应时间常数落在毫秒量级就够了比求解器的最小步长略大即可。第二个手段是把控制器单独放在一层。过热度 PI 控制器的输出是阀开度开度直接影响流进蒸发器的质量流量。如果控制器和物理模型混在一起解容易互相打架。我的做法是让控制器在每个求解器接受的时间步上更新一次用上一时刻的测量值相当于加了一个采样保持。这样物理方程和控制器解耦也顺便模拟了真实数字控制器的采样周期。注意虚拟容积法的容量系数如果取得过小会让雅可比矩阵的条件数急剧恶化取得过大又会把真实的压力动态掩盖掉。我一般从对应时间常数 1 毫秒起步然后看压力曲线的响应是否还能分辨出真实的物理特征。3. 换热器移动边界法和有限体积法的取舍3.1 移动边界法的方程推导与段长约束移动边界法的核心思路是换热器内部存在几个物理性质截然不同的区域过热、两相、过冷每个区域内部假设集中参数区域之间的边界位置随时间移动。这样你用很少的状态量就能抓住过热度在变、两相区在缩短这类关键信息。以蒸发器为例分两段两相段长度 L1和过热段长度 L2满足 L1 L2 L_total。每一段写质量和能量守恒$$ \frac{d}{dt}(\rho_i A_c L_i) \dot m_{in,i} - \dot m_{out,i} $$$$ \frac{d}{dt}(\rho_i A_c L_i u_i) \dot m_{in,i} h_{in,i} - \dot m_{out,i} h_{out,i} Q_i $$再给每段一个壁面能量方程$$ m_{w,i} c_w \frac{dT_{w,i}}{dt} Q_{ref,i} - Q_{air,i} $$这里藏着一个很容易被忽略的技术点段内的平均密度和平均内能怎么取。新手常犯的错是直接取进出口算术平均。在两相区这么干误差极大因为沿程干度是非线性增长的密度从饱和液到饱和气可能差两个数量级。正确做法是按空泡率沿程积分求平均或者至少用考虑滑移的平均空泡率公式。我实测过同一工况下算术平均和积分平均算出来的两相段质量能差 30% 以上直接导致压力响应速度完全不对。def evaporator_rhs(t, y, params, props): p, h1_out, h2_out, Tw1, Tw2, L1, L2 unpack(y) A, D, L_tot params.A, params.D, params.L_total # 段内平均物性按空泡率积分不能用进出口算术平均 rho1 props.two_phase_mean_rho(p, h_in, h1_out) u1 props.two_phase_mean_u(p, h_in, h1_out) m1 rho1 * A * L1 dm1 m_in - m_mid # 两相段质量守恒 dU1 m_in * h_in - m_mid * h1_out Q1 # 两相段能量守恒 # 由 dm1、dU1 反解 dp/dt 与 dh1_out/dt链式法则 2x2 线性方程组 J props.jacobian_mass_energy(p, h1_out) dp_dt, dh1_dt solve_2x2(J, np.array([dm1, dU1])) dL1 ... # 由过热段质量守恒反解 dL2 -dL1 # 总长度约束 return np.array([dp_dt, dh1_dt, dh2_dt, dTw1, dTw2, dL1, dL2])上面这段是伪代码级别的示意真实的实现里solve_2x2那个地方是移动边界法的心脏——你必须有能力把质量变化率和能量变化率转换成压力变化率和出口焓变化率。这需要用到状态方程对压力和比焓的偏导数也就是常说的物性导数。这部分如果拿不到解析导数靠有限差分去凑精度和速度都会很难看后面第 5 节会细说。3.2 段长趋零时的模式切换比方程本身更难写移动边界法真正折磨人的地方不在方程推导而在段长跑到零的那一刻。蒸发器过热度趋近于零的时候过热段长度 L2 会缩到 0然后理论上应该退化成单段全是两相模型状态量少一个。再过一会儿如果工况反过来过热段又得重新长出来。这一段逻辑如果你处理不好仿真会出现三种典型症状一是段长变成负数然后物性查询报错二是段长在 1e-5 附近反复抖动求解器步长被压到极小看起来像卡死三是切换之后能量不守恒曲线出现跳变。我的处理方式是三步走。第一步在状态方程里给段长做下限截断但截断不是硬来而是用一个光滑的过渡函数把趋近零的过程抹平def smooth_floor(x, eps1e-4, width1e-5): 平滑下限x 远大于 eps 时返回 x接近 eps 时平滑过渡 return 0.5 * (x eps np.sqrt((x - eps)**2 width**2))第二步用事件函数监视段长当它真的穿过阈值时终止当前积分由外层循环重建状态向量、切换模型模式、重新初始化再接着积分。这里必须用终止型事件不能让求解器带着已经失效的模型继续走。def event_superheat_vanishes(t, y, params, props): L2 y[-1] return L2 - 1e-4 event_superheat_vanishes.terminal True event_superheat_vanishes.direction -1 # 只在下降穿过时触发第三步切换之后做一次守恒量重建新的状态向量要保证总质量、总能量和切换前一致不能简单地把少的那个状态补零了事。这一步做不好你会看到压力在切换瞬间跳一下很多人以为是数值问题其实是建模问题。3.3 有限体积法分段数与数值扩散的平衡移动边界法不是万能的。它的假设是每一段内部充分混合、集中参数这在段比较长、沿程梯度不大的时候没问题一旦遇到除霜、或者制冷剂流量极低导致沿程干度变化剧烈集中参数的误差就藏不住了。这时候就需要有限体积法作为备选。有限体积法的做法是把换热器沿程切成 N 个控制体每个控制体写质量和能量守恒状态是每个控制体的比焓加一个全长度统一的压力。它的好处非常明显不需要模式切换。过冷、两相、过热这些区域的边界会自己浮现出来你不用写任何切换逻辑代码量反而更少。代价是状态数。一个换热器切 20 段冷凝器就是 20 个焓加 1 个压力再加 20 个壁温41 个状态两个换热器加起来八九十加上其他部件就上百了。这个规模对隐式求解器不算什么但对你调试来说就很痛苦——出了问题你根本不知道该看哪个状态。还有一个隐性的坑是数值扩散。一阶迎风格式在分段数不够的时候会把两相区和过热区之间的边界抹得很宽表现出来就是过热度好像一直在缓慢漂移找不到明确的拐点。我试过同一个工况下 10 段和 40 段的差别出口温度响应时间能差 15%。分段数不是越多越好但少于 15 段基本别指望能抓住区域边界的移动。3.4 两相区的空泡率与换热关联式选择顺序有讲究两相区的建模有两个必须做选择的地方空泡率怎么算换热系数怎么算。这两个选择互相独立但都会显著影响结果。空泡率方面均相模型假设气液同速算起来最省事但在低质量流速下会明显低估气相所占体积。滑移模型考虑了气相比液相跑得快更接近实际代价是需要一个滑移比的关联式。我的建议是如果你主要关心压力和过热度的宏观响应均相模型够用如果你要算两相段的质量存量、或者关心低流量下的动态就必须上滑移模型否则两相段质量会被系统性低估压力响应算出来偏快。换热系数方面单相区用经典的 Re^0.8 Pr^n 形式关联式就行成熟且稳定。两相区要区分沸腾和凝结因为两者机理不同蒸发器里是核态沸腾加对流蒸发冷凝器里是膜状凝结。这两类关联式在文献里有一大堆选哪个没有绝对答案但有两条实操经验一是优先选形式简单、导数连续的复杂关联式里常有分段和指数项容易把雅可比搞出尖角二是做完之后一定要用总换热量和设计值对一遍如果偏差超过 15%先怀疑关联式别急着调别的参数。选择项均相模型滑移模型建议两相段质量存量偏低 10%~30%接近实测关心压力动态必选后者计算成本低中差别不大导数连续性好一般需要额外光滑化适用流量范围高质量流速全范围低流量场景慎用均相3.5 我在两类场景下的实测取舍跑了一段时间之后我的结论是别只写一种两种都留着按场景选。做控制策略验证、启停过程、变负荷响应这一类用移动边界法。状态少、物理清晰、曲线干净你能一眼看出过热度是怎么被阀开度拽回来的。这类场景里段长很少真的跑到零切换逻辑大部分时间在待命。做除霜、做极端工况比如流量掉到设计值 10%、做多回路并联用有限体积法。这些场景里区域边界会大幅度移动甚至消失移动边界法的切换逻辑会被触发得非常频繁每次切换都带来一次重构误差累积起来很可观。有限体积法虽然状态多但省掉了这部分麻烦。4. 压缩机与节流让求解器不炸的写法4.1 效率查表与插值面构造压缩机在框架里是准稳态部件——它的机械和电气时间常数远小于热工过程所以不需要写微分方程直接用代数关系给定吸气状态和转速输出排气状态、质量流量、功率。核心是两个效率容积效率和等熵效率。厂家的样本上一般给的是不同蒸发温度、冷凝温度下的性能表。我见过有人直接拿样本上的制冷量反算流量这在小范围内可以但一旦工况超出样本范围外推出来的结果会很离谱。更稳的做法是用样本数据构造两个效率的二维插值面自变量是压比和转速然后在仿真里通过效率反算流量和功率。from scipy.interpolate import RectBivariateSpline class CompressorMap: def __init__(self, pr_grid, rpm_grid, eta_v_table, eta_s_table): # 输入必须是按压比排序、单调递增的网格 self.eta_v RectBivariateSpline(pr_grid, rpm_grid, eta_v_table, kx3, ky3) self.eta_s RectBivariateSpline(pr_grid, rpm_grid, eta_s_table, kx3, ky3) def mass_flow(self, p_suction, p_discharge, rpm, props): rho_s props.rho(p_suction, self.h_suction) pr p_discharge / p_suction eta_v float(self.eta_v(pr, rpm)) V_d self.displacement * rpm / 60.0 return eta_v * rho_s * V_d这里有几个容易翻车的细节。第一三次样条在网格边界附近容易出现非物理的振荡效率可能算出来大于 1 或者负数。要么降阶用线性插值要么在边界外面套一层裁剪。第二压比在启动瞬间可能小于 1还没建立压差这时候插值会跑到网格外必须有兜底逻辑比如强制压比下限为 1.0。第三效率查表本身不连续的话雅可比会出问题所以输入网格最好是单调的、密度足够的。4.2 节流孔口方程和过热度闭环节流装置有两种常见形式毛细管和电子膨胀阀。毛细管的特点是长度固定、流量由进出口压差和内部两相流动决定建模要写沿程摩擦压降比较麻烦电子膨胀阀简单得多用孔口方程就够$$ \dot m C_d A_{valve}(u) \sqrt{2 \rho_{in} (p_{in} - p_{out})} $$其中开度面积 $A_{valve}(u)$ 是阀门开度 u 的函数通常用厂家给的流量特性曲线或者简化为线性的流通面积关系。注意根号里如果是负数压差反了必须钳到零附近否则会得到复数。def valve_mass_flow(u, p_in, p_out, rho_in, Cd0.62): dp max(p_in - p_out, 0.0) A A_max * valve_curve(u) # valve_curve 在 [0,1] 内单调 return Cd * A * np.sqrt(2.0 * rho_in * max(dp, 1e-6))过热度闭环是这套框架里最能体现瞬态价值的部分。真实的电子膨胀阀控制器是一个 PI或者 PI 加点前馈控制器输入是蒸发器出口过热度和目标过热度的差输出是开度。这个闭环一旦接上整个系统的动态特性就变了——蒸发压力的响应不再只由压缩机决定而是被阀的动作牵着走。实现上要注意两件事。一是积分饱和开度被限制在 0 到 1 之间如果积分项一直累积退出饱和时会有一段严重的超调。加上抗饱和逻辑简单的做法是当输出饱和时停止积分。二是采样周期真实控制器有几十毫秒到几百毫秒的采样周期如果你让控制器每个求解器内部步都更新就等于给了它一个无限快的响应算出来的过热度会过于平滑。我一般用一个固定的控制周期作为采样保持跟真实设备对得上。4.3 硬切换改成软过渡三种我用过的手段制冷系统里到处是开关压缩机启停、风机启停、电磁阀开闭、除霜切换。这些东西在数学上是不连续点处理不好就会让求解器在那个时刻反复尝试、步长越压越小。第一种手段是把阶跃改成斜坡。压缩机从 0 到 3000 rpm不要在一瞬间完成给一个 0.5 到 2 秒的斜坡。物理上这也不失真——真实变频器也有加速时间。这一步能解决掉大概一半的切换问题。第二种手段是用光滑的 min/max 替代硬判断。比如热流方向、压差的正负号用硬 if-else 写会得到不可导的右端项。换成平滑版本def smooth_max(x, y, k1e3): return 0.5 * (x y np.sqrt((x - y)**2 1.0/k**2))参数 k 要选到物理上该切换的地方确实切换了但过渡区内不是极端陡峭。k 太大等于没抹平太小会把有意义的物理突变抹没了。第三种手段是用事件 状态重置而不是让方程自己去穿过不连续点。凡是离散状态改变的事情比如压缩机停机、风机切档都应该做成终止型事件停下来改模型再重新开始积分。这样做的好处是每一步都是光滑的求解器不会在一个尖叫点上挣扎。5. 物性调用最容易把仿真拖垮的地方5.1 反查迭代不收敛的四种成因用压力加比焓做状态意味着你在每一次右端项计算时都要反查一堆物性密度、温度、熵、内能、饱和线上的干度。这些反查在单相区是直接查表插值很快在两相区和临界区附近就是迭代求解慢而且可能不收敛。我踩过的收敛问题基本能归成四类。第一类是超出了有效范围。比如压力算出来是个负数或者超过临界压力你还在按两相区去查饱和线必然失败。第二类是输入本身是错的压力正常但比焓落在液区深处你却按两相区逻辑去算干度出来的干度是负的再拿负数去算空泡率就崩了。第三类是物性库的默认迭代容差不够在临界点附近残差下不去的问题被放大。第四类是求解器的试探步给了非物理的状态。隐式方法在迭代中间会尝试各种离谱的值这些值本身没有物理意义但物性函数不能崩崩了求解器就完蛋。对应的处理办法很朴素所有物性调用都要有边界保护和兜底返回。压力钳到 [p_min, p_crit0.999]温度钳到 [T_min, T_crit0.999]干度钳到 [0, 1]。听上去很脏但这是工业级仿真的常规做法。你在迭代中间给一个外推值比让整个求解崩掉要好得多。当然用完之后要记录有多少次触及了边界这个计数能帮你判断仿真是不是在硬撑。5.2 导数不连续是怎么把隐式求解器毁掉的隐式求解器每一步都要解一个非线性方程组靠牛顿迭代需要的核心是雅可比矩阵。如果物性函数在某处不可导或者导数跳变牛顿迭代就会失去二次收敛性表现出来是步长被压缩到极小、迭代次数暴涨、最后报步长太小或者收敛失败。最典型的不可导点就是饱和线。你在两相区算密度沿着等压线走干度到 1 的时候流体刚好变成饱和气密度的斜率有一个转折。如果你的状态变量刚好在饱和线附近来回穿越求解器每一步都在踩这个转折。过热度在 0 附近的情况是最常见的现场因为蒸发器出口本来就该在饱和气或微过热状态。处理这个问题的办法有几个层次。最省事的是加滞环定义过热度时用一个小的死区让状态在 ±0.2 K 以内不改变判定逻辑。这不能消除不可导但能避免反复穿越。更彻底的是光滑化把分段的物性用一个连续的过渡函数接起来代价是引入一点非物理的平滑。最讲究的是把相区作为离散模式用事件切换切换之后重新初始化这样每段积分内部都是光滑的。第三种工作量最大但在要求高的场景里是唯一可靠的做法。提示判断仿真是不是踩了不可导点一个很实用的小技巧是记录每个接受步的雅可比条件数和牛顿迭代次数。如果这两个数在某个时间段突然集体变大八成就是某个状态穿过了饱和线或者区域边界。5.3 缓存、表格与导数提供物性调用是这类仿真的第一性能瓶颈没有之一。我统计过在没做任何优化的情况下一次完整的启停仿真里物性函数的调用次数占整个运行时间的 60% 到 80%。优化路径有三条收益从大到小排列。第一条是缓存。同一个 (p, h) 组合在雅可比数值差分的时候会被反复查询——差分本身就要在附近扰动好几次。做一个带 LRU 的缓存命中率通常能到 40% 以上直接省掉这部分时间。from functools import lru_cache lru_cache(maxsize200_000) def _rho_ph(p_key, h_key): # p_key, h_key 是四舍五入后的整数保证缓存能命中 return _backend_rho(p_key * 1.0, h_key * 1.0) def rho_ph(p, h): return _rho_ph(round(p, 1), round(h, 1))注意这里的取整精度要拿捏取整太粗物性精度不够太细缓存命中率低。我一般把压力取到 0.1 Pa、比焓取到 0.1 J/kg对工程精度完全够用。第二条是预建表格 双线性插值。把整个工作范围压力从 0.1 到 5 MPa比焓从液区到气区预先离散成一张二维表运行时直接插值。这一步能把单次查询从微秒级降到纳秒级提速一两个数量级。代价是内存和精度表太粗会丢掉饱和线附近的细节。我的做法是在饱和线附近加密网格远离饱和线的地方用粗网格。第三条是提供解析导数。如果你的物性后端能给出对压力的偏导、对比焓的偏导那移动边界法里那个 2x2 的线性方程组就不用数值差分了精度和速度双赢。这条的收益在状态数多的时候才明显状态少的时候不如前两条划算。6. 求解策略刚性、事件与一致性初始化6.1 为什么显式方法一步都走不动前面提过刚性比在 10^3 到 10^4这里给一个具体的量感。假设最快的时间常数是管壁的 1 秒最慢的是房间的 3000 秒。显式欧拉法的稳定步长要求大约是 2 倍的最快时间常数也就是 2 秒。听起来还行问题是这只是稳定性的必要条件精度上你还得再小一个量级。而且一旦你的模型里加进了压力节点的虚拟容积时间常数毫秒级显式步长立刻掉到毫秒以下算一小时的工况就是 360 万步。隐式方法BDF、Radau 这类没有稳定性对步长的硬约束它可以根据精度自适应地放大步长。同一个工况隐式方法可能几千步就走完了。这就是为什么框架从第一天起就必须按隐式设计——不是性能优化是可行性问题。from scipy.integrate import solve_ivp sol solve_ivp( funrhs, t_span(t0, t1), y0y0, methodBDF, # 刚性系统首选Radau 也可以但每步更贵 rtol1e-5, atol1e-7, events[ev_superheat_vanishes, ev_pressure_limit], max_stepnp.inf, # 交给自适应别自己限死 dense_outputFalse, )参数上有两个经验rtol不要一味收紧1e-5 到 1e-6 对工程问题足够了收到 1e-9 只会让步长变小、时间变长曲线看不出区别。max_step千万别随手设一个很小的值那等于把自适应步长的优势废掉了只在需要强制输出采样点的时候才用。6.2 事件函数与状态重置事件是这套框架里最需要小心设计的部分。事件的本质是在某个条件成立的那一刻停下来但停下来的原因不同后续处理完全不同。我把事件分成三类分别处理事件类型例子处理方式物理边界压力超过安全阈值、段长趋零终止积分切换模式或报错控制动作压缩机启停、阀开度阶跃终止积分修改输入继续输出采样每 1 秒记录一次曲线非终止仅记录写事件函数的时候有两个坑。一是符号方向。事件函数写成L2 - threshold如果只关心它从正变负就要设direction-1否则求解器会把它在正向穿越时也当成事件导致莫名其妙的提前终止。二是事件的密集触发。如果事件函数在零附近抖动会连续触发终止外层循环陷入死循环。防抖的办法是给事件加一个最小间隔或者触发后强制推进一小段时间再重新布防。状态重置是另一处高危区。切换模式之后新的状态向量维度和含义都变了如果你只是简单地补零或者截断质量、能量就不守恒了压力会出现一个物理上不存在的跳变。正确的做法是从旧状态反算出各段的存量、总能量然后构造新状态使得这些总量保持不变。这一步的代码不长但必须写。6.3 从稳态解起步的一致性初始化瞬态仿真最常见的启动失败原因是初始条件不自洽。你随手给一个压力、一个焓它们在物理上可能对应一个完全不平衡的状态求解器一上来就要处理巨大的残差直接发散。正确的做法是先解稳态再拿稳态解做初值。具体来说固定压缩机转速、阀开度、外部工况把所有的 dy/dt 置零解一个非线性方程组。这个方程组的规模跟状态数一样用牛顿法解很快。def solve_steady_state(y_guess, params, props, tol1e-8, max_iter50): y y_guess.copy() for _ in range(max_iter): F steady_residuals(y, params, props) # 所有 dy/dt 0 if np.linalg.norm(F, np.inf) tol: return y, True J numerical_jacobian(steady_residuals, y, params, props) dy np.linalg.solve(J, -F) # 阻尼牛顿限制单步变化量避免跳到非物理区域 scale min(1.0, 0.5 / max(np.max(np.abs(dy)), 1e-12)) y y scale * dy return y, False这段里有几个细节很关键。阻尼不能省纯牛顿在物性高度非线性的系统里几乎必然过冲。初值的选择也重要最好从设计工况附近起步别从一个完全不相干的点开始。收敛判据用无穷范数而不是二范数因为你要保证每个方程都满足而不是总体上凑合。稳态解出来之后还有一步检查它是否物理合理。看四个数——蒸发压力、冷凝压力、过热度、过冷度跟设计值对一遍。如果对不上别急着往下跑瞬态先把稳态模型调对。这一步能省掉后面几天的排查时间。6.4 雅可比的稀疏结构与解析导数状态数上了几十以后雅可比矩阵的数值差分成本就上来了。默认的有限差分要对每个状态各扰动一次一次雅可比就是 N 次右端项计算而右端项里最贵的是物性查询。N50 的时候一次雅可比可能就要几十毫秒几步下来仿真时间就爆了。两个优化方向。一是利用稀疏结构。部件之间的耦合是通过端口发生的一个蒸发器的状态只跟它相邻的几个部件有关雅可比天然是块稀疏的。把这个稀疏模式告诉求解器能省掉大量零元素的差分。二是提供解析或半解析导数。物性那一层的偏导如果后端能给就用后端的给不了至少把线性部分的导数写出来只对非线性最强的物性部分用差分。我实测过的对比一个 60 状态的模型纯数值差分跑完 600 秒仿真要 8 分钟加上稀疏模式和物性缓存之后降到 50 秒左右。这个差距在参数扫描的时候会被放大几十倍——你要扫 100 组工况8 分钟就变成 13 个小时。7. 怎么判断曲线是对的验证与一次真实排错7.1 先对齐稳态设计点再谈瞬态瞬态模型最尴尬的地方是它没有直接的可比对象。稳态有选型软件可以对照瞬态你只有一条自己算出来的曲线。所以验证必须分层做。第一层是稳态退化验证。把所有导数项置零模型应该退化成一个稳态循环模型。这时候算出来的换热量、压比、过热度、过冷度应该跟独立的选型计算结果对得上。我给自己定的容差是换热量和功率 ±5%压力和温度 ±2%。达不到这个说明基础模型有问题别往瞬态上走。第二层是能量守恒自检。整个循环里压缩机的输入功应该等于冷凝器放热减去蒸发器吸热忽略向环境的散热损失。这个残差在每一个时间点都应该接近零。def energy_residual(sol): W_comp sol[compressor_power] Q_cond sol[condenser_heat_rejection] Q_evap sol[evaporator_heat_absorption] return (W_comp Q_evap - Q_cond) / max(Q_cond, 1e-6)这个残差如果在 0.5% 以内说明你的能量方程写得没大问题。如果超过 2%一定是某个地方的正负号错了或者某个流量项漏了。第三层才是动态特性对比。如果有实验数据比对阶跃响应的时间常数、超调量、稳态误差。没有实验数据的话至少做两个合理性检查响应方向对不对转速升高排气压力应该升高时间常数在不在物理直觉的范围内管壁主导的响应不该是毫秒级。7.2 一次蒸发器出口温度锯齿的完整排查链路讲一个我实际遇到的、很有代表性的问题。现象是仿真跑到 200 秒左右蒸发器出口温度开始出现高频锯齿幅度不大0.3 K 左右但明显不是物理现象而且锯齿出现之后求解器步长明显变小仿真变慢。我的排查顺序是这样的。第一步先排除求解器容差。把 rtol 从 1e-5 收到 1e-8重新跑。锯齿幅度没有变化只是变密了一点。这一步排除了精度不够的可能。第二步检查是不是控制器的采样周期问题。把 PI 控制器的采样周期从 100 ms 改成 10 ms锯齿频率跟着变了。这说明锯齿跟控制器有耦合但根因不一定在控制器。第三步加日志记录每次物性查询的输入输出。在饱和线附近加一个探针记录每次查询的 (p, h) 和返回的干度。跑完之后一看问题清楚了出口状态在 h 略大于饱和气焓的位置反复穿越干度在 0.999 和 1.001 之间跳。当干度小于 1 时走两相区物性路径大于 1 时走单相区路径两条路径在衔接处的密度算出来有 0.5% 的不连续。第四步定位到物性路径切换。这就是典型的饱和线不可导。控制器在过热度接近目标值的时候会把阀开度微调导致出口状态刚好在饱和线上下浮动每次浮动都触发一次物性路径切换。第五步修复。三个改动叠加上去一是物性查询在两相/单相边界处加一个过渡带宽度 1e-3 的干度带内用加权平均二是过热度计算加 0.1 K 的滞环避免控制器在目标值附近过度动作三是控制器的输出加一阶滤波时间常数 0.2 秒对应真实阀门的机械响应。第六步验证。锯齿消失能量残差保持在 0.2% 以内出口温度曲线变得平滑同时 200 秒附近的真实动态特征过热度从 6 K 回落到 5 K 的过程仍然保留。求解器步长恢复到正常水平。这个案例给我的教训是高频小幅振荡八成是某个开关在反复动作。物性路径切换、相区判定、控制器饱和、事件触发这几个是重点怀疑对象。思路就是一层层剥先排除求解器再排除控制器最后落到模型本身。7.3 几个我用来快速判断哪里错了的启发式跑得多了之后会形成一些直觉这里分享几条我觉得最有用的。压力曲线不合理地平滑——通常是虚拟容积的容量系数取太大了把真实动态掩盖了。把系数调小一个量级试试。压力曲线有台阶状——通常是某个代数环没解开压力在迭代中被卡在某个值上。检查是不是有部件在算流量时用了上一时刻的压力。过热度一直贴着 0——不是模型错是阀开度给太大了或者压缩机排量相对换热器太小。先检查设计点。能量残差随时间缓慢漂移——通常是有个地方用了显式更新比如壁温用上一时刻的值算了热流但没进能量方程。查一下每个能量方程里的热流项是不是都用当前状态算的。仿真在某一个固定时刻附近总是失败——看那个时刻有没有事件触发。事件的符号方向或者防抖没写对。8. 从能跑到好用工程化与加速的优先级8.1 组件注册与配置驱动第一版能跑通之后最想做的是把各种参数硬编码清理掉。我采用的方式是组件注册加配置驱动每个部件类注册一个名字配置文件里描述回路拓扑和参数运行脚本只负责读配置、建对象、连端口、调求解器。REGISTRY {} def register(cls): REGISTRY[cls.__name__] cls return cls register class EvaporatorMB(Component): ... def build_system(cfg): comps {name: REGISTRY[c[type]](name, **c[params]) for name, c in cfg[components].items()} for a, b in cfg[connections]: connect(comps[a].outlet, comps[b].inlet) return comps这么做的好处在做参数扫描的时候体现得最明显。以前换个压缩机要改代码现在改一行配置就行以前跑 100 组工况要写循环现在写个脚本生成 100 份配置扔到多进程池里就行。这一步的投入产出比非常高我建议在第一版跑通之后立刻做别等到代码乱了再重构。8.2 加速手段的收益排序别把力气花错地方性能优化最容易犯的错是先去做那些看起来很酷但收益很小的事。按我的实测收益从大到小是这样的手段典型加速比实现成本建议时机物性缓存2~5 倍低立刻做物性表格化 插值10~50 倍中状态数 30 时做稀疏雅可比 解析导数3~10 倍中高单次仿真 1 分钟时做减小状态数降阶、合并5~20 倍高有明确精度验证后做JIT 编译热点函数2~5 倍中前三项做完还不够时多进程参数扫描近似线性低需要扫工况时立刻做注意最后一条和前面几条不是一个维度前面是让单次仿真变快最后一条是让一批仿真变快。如果你的使用场景是参数扫描或者优化多进程的收益往往比任何单次优化都大而且实现成本最低几行代码的事。还有一个常被忽略的点输出频率。默认的求解器只在内部步长上输出如果你设了很密的输出点求解器会被迫缩小步长去对齐白白浪费。正确做法是让求解器自由步进用一个插值或者独立的采样循环去生成等间隔输出。8.3 后面还能往哪加以及我为什么暂时不加框架跑通之后想加的东西会排成长队多回路、混合工质、润滑油循环、结霜除霜、风扇动态、变频器损耗、噪声。我的做法是每个都问一句它会不会改变我要验证的那个结论。结霜会影响蒸发器换热和风量如果你做的是冬季热泵长时间运行这个必须加如果只做启停几秒钟的动态可以暂时忽略。风扇动态的时间常数通常在秒级跟管壁接近加进来会让模型更真实但也会让状态数翻倍。润滑油主要影响低流量下的回油和换热跟压力动态关系不大。我的判断标准一直是先把一个闭环做到自己敢信再往上面叠功能。一个能解释清楚每一条曲线为什么长这样的简单模型比一个功能齐全但没人说得清的黑箱有价值得多。这套框架后面怎么长取决于你拿它解决什么问题而不是取决于它能支持多少特性。最后分享一个我自己用得很顺的小习惯每次改动模型之后不要只看新算出来的曲线而是把改动前后的两条曲线画在同一张图上同时把能量残差和求解器步长也画出来。曲线告诉你结果变了残差告诉你守恒有没有破步长告诉你数值上是不是变难了。这三个信息放在一起绝大多数改动的对错一眼就能判断比反复读代码快得多。