pipe-system-simulation-test/docs/superpowers/specs/2026-04-14-0d-1d-tank-pipe-blowdown-mvp-design.md

0D–1D 气瓶–管路耦合瞬态放气仿真 MVP

日期：2026-04-14 状态：草案，待评审 作者：He Yunqin（通过 Claude Code / superpowers:brainstorming）

---

1. 目标与范围

1.1 目标

实现一个可运行的瞬态仿真程序，模拟以下场景：

• 高压气瓶 V₁ = 5 m³，P₁ = 10 MPa，T₁ = 300 K
• 低压气瓶 V₂ = 10 m³，P₂ = 0.101325 MPa（1 atm），T₂ = 300 K
• 连接管路 L = 1 m，D = 5 mm，划分为 N = 20 个有限体积单元
• 工质：理想气体，γ = 1.4，R = 287 J/(kg·K)
• 初始条件：管路与低压瓶同压同温（P₂, T₂），高压瓶独立
• 仿真时长：t_end = 0.1 s（捕获开启瞬间的激波与膨胀波过程）

程序应输出：

1. 双瓶 P(t)、T(t) 时间曲线图
2. 管内 P(x)、u(x)、T(x) 演化动画（GIF）
3. 完整时间序列数据文件（.npz），支持离线加载任意时刻任意空间点的状态

1.2 MVP 范围边界

本 MVP 做：

• 0D 集中参数气瓶 + 1D 可压缩欧拉方程管路
• 虚网格（Ghost Cell）+ HLL Riemann 求解器耦合
• 一阶空间重构 + 一阶显式欧拉时间推进
• CFL 动态步长控制
• 基本 pytest 单元测试（共 10 条）

本 MVP 明确不做：

• 摩擦源项（Darcy–Weisbach 或 Fanning）
• 壁面传热（绝热假设）
• 二阶空间格式（MUSCL / limiter）
• 二阶时间格式（SSP-RK2/RK3）
• 局部阻力修正（入口/出口损失系数 ζ）
• HLLC、Roe、exact Riemann 等其他格式
• 网格收敛性扫描
• 真实气体状态方程
• 多管网络、分叉、汇合

---

2. 物理与数值模型

2.1 0D 气瓶模型

对每个气瓶，假设内部状态均匀、气体静止，采用质量与总内能两个守恒量：

dm/dt = ṁ_in                                      (质量守恒)
dU/dt = Ḣ_in  (= ṁ·h_t,in)                        (能量守恒，开口系第一定律)

其中：

• m = 总质量 [kg]
• U = 总内能 [J]，对理想气体 U = P·V/(γ−1) = m·R·T/(γ−1)
• Ḣ = 总焓流 [W]，正值表示流入

派生量通过状态方程按需计算：

ρ = m/V
T = (U/m) · (γ−1)/R
P = ρ·R·T = U·(γ−1)/V

设计取舍：选 (m, U) 作底层状态而非 (P, T) 的理由见 §4.3。

2.2 1D 管路模型

一维可压缩欧拉方程（守恒形式）：

∂W/∂t + ∂F(W)/∂x = 0

W = [ρ, ρu, ρE]ᵀ
F = [ρu, ρu² + P, u·(ρE + P)]ᵀ

其中 E = e + u²/2，e = P/(ρ(γ−1))。

离散化：有限体积法，cell-averaged piecewise constant（一阶重构）。

W_i^{n+1} = W_i^n − (Δt/Δx)·(F_{i+1/2} − F_{i−1/2})

• N = 20 个单元
• Δx = L/N = 0.05 m
• 单元中心坐标 x_i = (i + 0.5)·Δx, i = 0..19

2.3 耦合：虚网格 + HLL

虚网格定义：管路左右各增加一个"虚拟单元"，每个时间步开头根据当前气瓶状态填充：

W_ghost_L = [ρ₁, 0, P₁/(γ−1)]ᵀ                   ← 高压瓶
W_ghost_R = [ρ₂, 0, P₂/(γ−1)]ᵀ                   ← 低压瓶

速度取 0 的假设：把气瓶视为滞止状态（stagnation state）的无穷大储气罐，真实的加速过程由 HLL 在界面上解出。这是建模简化，不是物理真相；其代价是瓶内到界面的速度-压力关系并非严格等熵，对短管、大压差场景误差可接受。

注意区分两件事：

1. 系统总能守恒：从"通量双用"的数学结构直接推出（见 §2.4），与 HLL 精度无关，严格到机器精度成立。
2. 气瓶能量更新的物理一致性：HLL 的 flux[2] = u(ρE + P) 等同于总焓流密度 ρu·h_t。对于 u_ghost=0 的滞止态气瓶，h_t,ghost = c_p·T_tank，这恰好是开口系能量方程对气瓶的正确源项——即便 Riemann 求解器引入数值耗散，能量守恒量也被 HLL 的守恒形式严格守住。

HLL 数值通量（Harten–Lax–van Leer）：

S_L = min(u_L − a_L, u_R − a_R)
S_R = max(u_L + a_L, u_R + a_R)

          ┌ F_L                                          if S_L ≥ 0
F_HLL =   │ (S_R·F_L − S_L·F_R + S_L·S_R·(W_R − W_L))    if S_L < 0 < S_R
          │  ────────────────────────────────────
          │                S_R − S_L
          └ F_R                                          if S_R ≤ 0

其中 a = √(γP/ρ)。

关键性质：HLL 自动处理双向流、激波、膨胀波和截流；无需手动判断流动方向或是否达到声速。

2.4 守恒性与通量双用

核心耦合原理：每个时间步，HLL 在管路左右两个边界各给出一个数值通量向量。同一组通量被两边共享：

• 用于更新管路首/末单元的保守变量（有限体积更新）
• 乘以管路截面积 A 后，作为气瓶的 ṁ, Ḣ 源项更新气瓶

这一"通量双用"机制从数学上天然保证系统总质量和总能量守恒到机器精度，是验证实现正确性的唯一充分必要条件。

---

3. 时间推进算法

3.1 单步推进的 7 个阶段

阶段 1 — 计算 CFL 时间步长
  a_i   = √(γ P_i / ρ_i)                 for each cell i
  a_max = max over i of (|u_i| + a_i)
  Δt    = CFL · Δx / a_max
  Δt    = min(Δt, t_end − t)             夹到 t_end

阶段 2 — 冻结气瓶状态为虚网格
  W_ghost_L = tank1.ghost_state()       # 本步用的高压瓶快照
  W_ghost_R = tank2.ghost_state()       # 本步用的低压瓶快照
  （管路 W 的快照留给 pipe.step 内部处理，solver 不碰）

阶段 3 — 计算两个边界通量（solver 层）
  flux_L = HLL(W_ghost_L,     pipe.W[:,0])              左边界 (W_ghost_L ↔ cell 0)
  flux_R = HLL(pipe.W[:,N-1], W_ghost_R)                右边界 (cell N-1 ↔ W_ghost_R)

  注：只有这两个"边界通量"由 solver 计算，因为它们需要知道 ghost state。
      所有 N-1 个"内部通量"由 pipe.step() 内部计算，solver 不参与。

阶段 4 — 管路一步推进（pipe.step 内部完成）
  pipe.step(flux_L, flux_R, Δt):
      W_snap   = pipe.W.copy()
      flux_int = [HLL(W_snap[:,i-1], W_snap[:,i]) for i in 1..N-1]    内部 N-1 个通量

      首单元:   pipe.W[:,0]   = W_snap[:,0]   − (Δt/Δx)·(flux_int[1]   − flux_L)
      末单元:   pipe.W[:,N-1] = W_snap[:,N-1] − (Δt/Δx)·(flux_R        − flux_int[N-1])
      中间:     pipe.W[:,i]   = W_snap[:,i]   − (Δt/Δx)·(flux_int[i+1] − flux_int[i])
                                                              for i in 1..N-2

阶段 5 — 同一组边界通量更新气瓶
  fL_A = flux_L · A                          # (3,) 向量 × 标量面积
  fR_A = flux_R · A
  tank1.apply_flux(mdot=fL_A[0], edot=fL_A[2], dt=Δt, sign=-1)    流出
  tank2.apply_flux(mdot=fR_A[0], edot=fR_A[2], dt=Δt, sign=+1)    流入

阶段 6 — 推进时间
  t += Δt

阶段 7 — 记录历史
  history['t'].append(t)
  history['P1'].append(tank1.P)
  history['T1'].append(tank1.T)
  history['P2'].append(tank2.P)
  history['T2'].append(tank2.T)
  history['W_hist'].append(pipe.W.copy())

3.2 关键顺序约束

1. 必须先一次性算完所有界面通量，再一次性更新所有单元。如果边算边更新，会导致"新 W 和旧 W 混用"。
2. 管路和气瓶的更新使用同一组冻结通量。若用已被更新的 tank 状态重算左边界通量，即违反守恒性。
3. 阶段 4 和阶段 5 的先后顺序不重要（彼此独立），但不能交错。

3.3 CFL 估算与步数预估

初始时刻（管路全场 = 低压瓶态）：

a₀ = √(1.4 · 101325 / 1.177) ≈ 347 m/s
u₀ = 0
Δt₀ = 0.5 · 0.05 / 347 ≈ 7.2 × 10⁻⁵ s

激波进入管路后 |u| + a 可达 600–800 m/s，Δt 收缩到 ~30–40 μs。

t_end = 0.1 s 预计对应 约 2000–3000 步。NumPy 向量化 N=20 的数组，预计全程 <15 秒 跑完。

---

4. 组件与模块划分

4.1 文件布局

pipe-system-simulation-test/
├── src/
│   ├── config.py       # 常量与工况（γ, R, V, P_init, L, D, N, CFL, t_end）
│   ├── riemann.py      # HLL 数值通量（纯函数）
│   ├── tank.py         # Tank 类（0D，状态 = (mass, U)）
│   ├── pipe.py         # Pipe 类（1D 有限体积场 W, 推进方法）
│   ├── solver.py       # 时间循环与气瓶-管路耦合编排
│   ├── output.py       # save_history / plot_timeseries / make_animation
│   └── main.py         # 入口：装配、调度、sanity check
├── tests/
│   ├── test_riemann.py
│   ├── test_tank.py
│   ├── test_pipe.py
│   └── test_integration.py
├── conftest.py         # sys.path 注入，让 tests/ 能从 src/ 导入
├── results/            # 运行产物（.npz, .png, .gif）
├── cases/              # （保留，后续扩展多工况用）
├── docs/
│   └── superpowers/specs/2026-04-14-0d-1d-tank-pipe-blowdown-mvp-design.md  # 本文件
└── scripts/            # （保留，后续工具脚本用）

4.2 模块职责

模块	职责	对外接口
config.py	纯数据：物理常数、工况、仿真控制参数	顶层常量
riemann.py	HLL 数值通量计算	hll_flux(W_L, W_R, gamma) → ndarray(3,)
tank.py	0D 气瓶状态与演化	Tank 类，ghost_state(), apply_flux()，派生属性 P/T/rho
pipe.py	1D 管路状态与一步推进	Pipe 类，primitives(), max_wave_speed(), step()
solver.py	装配时间循环、调用 HLL、协调气瓶与管路更新	run(tank1, tank2, pipe, t_end, cfl, verbose=False, log_every=100) → history: dict
output.py	持久化与可视化	save_history(), plot_tank_timeseries(), make_pipe_animation()
main.py	入口脚本：构造对象 → 调 solver → 调 output → sanity check	main()

4.3 关键设计取舍

4.3.1 Tank 底层状态选 (mass, U) 而非 (P, T)

理由：

• 守恒律在 (m, U) 空间是线性的：m += ṁ·Δt，U += Ḣ·Δt，两次加法即完成，零公式展开
• (P, T, ρ) 三者由状态方程约束，只有 2 个自由度；若都当第一身份存，必须人工保证同步，任何一步漏更新就违反状态方程
• HLL 的 flux[2] = u(ρE + P) 恰好等于总焓流密度 ρu·h_t，乘以 A 后就是开口系能量方程的右端项 Ḣ；写成 dU/dt = Ḣ 直接对应代码 tank.U += flux[2]·A·dt，无需手工添加"流动功修正"
• 派生属性（tank.P, tank.T, tank.rho）通过 @property 实时计算，永远与 (m, U) 一致，外部读不到过时值

代价：每次访问 tank.P 需要一次除法 + 一次乘法。对 0D 气瓶完全可忽略。

4.3.2 界面通量"双用"，由 solver 显式编排

• solver.run() 在每步先算边界 flux，再把 flux 分发给 pipe.step(flux_L, flux_R, dt) 和 tank.apply_flux(ṁ, Ḣ, dt, sign)
• pipe.step() 内部负责所有内部界面通量的计算和管路单元更新
• 气瓶不知道管路存在，管路不知道气瓶存在；耦合知识集中在 solver 一处

4.3.3 ghost_state() 作为 Tank 的方法

• "如何把自身暴露成虚网格"是 tank 的职责，不是 solver 的职责
• solver 代码里就是 hll_flux(tank1.ghost_state(), pipe.W[:,0], γ)，可读性最高
• 未来如果要扩展成"带局部阻力修正的虚网格"，只改 Tank.ghost_state() 即可

4.3.4 sign 参数约定

Tank.apply_flux(mdot, edot, dt, sign) 接受一个显式的 sign 参数：

• sign = -1：气流"流出"瓶子。对 tank1（高压瓶），左边界通量向右为正，即流出，所以用 -1
• sign = +1：气流"流入"瓶子。对 tank2（低压瓶），右边界通量向右为正，即流入，所以用 +1

显式 sign 比在 Tank 内部判断"我是 tank1 还是 tank2"更干净。

---

5. 数据流与历史记录

5.1 History 数据结构

history = {
    't':      np.ndarray,     # shape (n_steps,)           时间序列
    'P1':     np.ndarray,     # shape (n_steps,)           高压瓶压力
    'T1':     np.ndarray,     # shape (n_steps,)           高压瓶温度
    'P2':     np.ndarray,     # shape (n_steps,)           低压瓶压力
    'T2':     np.ndarray,     # shape (n_steps,)           低压瓶温度
    'W_hist': np.ndarray,     # shape (n_steps, 3, N)      管路全程守恒变量
}

循环中用 Python list + append 收集，循环结束后一次性 np.stack / np.asarray 转成数组。

内存估算：

n_steps ≈ 3000, 3 vars × 20 cells × 8 bytes
→ 每步 480 B, 全程 ≈ 1.4 MB

完全无压力，每步都记，不做降采样。

5.2 持久化：results/history.npz

np.savez_compressed(
    "results/history.npz",
    t        = history['t'],
    P1       = history['P1'], T1 = history['T1'],
    P2       = history['P2'], T2 = history['T2'],
    W_hist   = history['W_hist'],
    x        = pipe_cell_centers,   # shape (N,)
    dx       = pipe.dx,
    gamma    = GAMMA,
    R_gas    = R_GAS,
    area     = pipe.area,
)

未来加载与访问示例：

d = np.load("results/history.npz")
rho = d['W_hist'][:, 0, :]                    # (n_steps, N)
u   = d['W_hist'][:, 1, :] / rho
P   = (d['W_hist'][:, 2, :] - 0.5*rho*u**2) * (d['gamma'] - 1)

5.3 可视化产物

results/
├── history.npz              全量时间序列，可离线复盘
├── tank_pressure.png        P1(t), P2(t) 双曲线
├── tank_temperature.png     T1(t), T2(t) 双曲线
└── pipe_animation.gif       管内 P(x), u(x), T(x) 三子图演化

动画细节：

• 用 matplotlib.animation.FuncAnimation + PillowWriter（GIF 输出，免 ffmpeg 依赖）
• 布局：3 个纵向子图，共享 x 轴（管路空间坐标 0–1 m）
• 每帧标题：t = X.XXXX s (step k/n_steps)
• 默认 stride = 10，约 300 帧，GIF 文件预计 5–10 MB
• 数据存储仍为全量；stride 只影响动画帧数

---

6. 错误处理

6.1 验证点清单

位置	检查	失败时
程序启动	results/ 目录存在（os.makedirs(exist_ok=True)）	静默创建
config.py 顶层	γ>1, R>0, V>0, L>0, D>0, N≥2, P>0, T>0, CFL∈(0,1]	AssertionError
时间步开始	管路所有单元 ρ>0, P>0	RuntimeError("负密度/负压力 at step k")
CFL 步长	Δt > 1e-12	RuntimeError("dt 退化")
HLL 内部	恢复原始变量时 ρ>0, P>0	ValueError + 打印 W_L, W_R
Tank 更新后	mass > 0	RuntimeError("气瓶质量非正")
仿真结束	n_steps > 0	RuntimeError("一步都没跑")
main.py 末尾	总质量相对误差 < 1e-10	AssertionError（作为 sanity check 打印）
main.py 末尾	总能量相对误差 < 1e-10	AssertionError（作为 sanity check 打印）

6.2 策略

• 异常直接往上抛，不 catch、不重试、不降级。让 Python traceback 定位问题
• 不 try/except 包装 main.py
• solver.run() 提供 verbose=False 开关（默认关）和 log_every=100 打印间隔，用于调试时查看步进日志

6.3 明确不做的检查

• 负质量流量判断（HLL 自动处理双向）
• 显式判截流（Riemann 自动处理）
• 动态 CFL 调整策略
• 并行/线程安全

---

7. 测试策略

7.1 测试清单（共 10 条，4 个文件，预计 <5 秒）

tests/test_riemann.py — HLL 求解器

1. 左右状态相同 → 返回纯物理通量（无数值耗散）
2. 静止接触间断（两侧 u=0，仅 ρ 不同）→ 质量通量和能量通量应为 0，动量通量 = P
3. Sod 激波管初值的符号与量级检查（不做精确值断言，检查 flux[0], flux[1], flux[2] 均 > 0）

tests/test_tank.py — 气瓶

4. 初始状态自洽：给定 (P, T, V) 构造 → mass 和 U 满足理想气体关系
5. ghost_state 是滞止态：动量项 = 0，能量项 = P/(γ−1)
6. apply_flux 方向性：sign=-1 应减质量、减内能、降压

tests/test_pipe.py — 管路

7. 均匀初始化：所有单元的 W 完全相同；primitives() 返回 u=0, P=P_init
8. 静止状态 + 一致压力通量 → 一步后 W 不变（机器精度）

tests/test_integration.py — 端到端守恒律

9. 短仿真总质量守恒：t_end=1e-3，|Δm_total/m_total| < 1e-10
10. 短仿真总能量守恒：同上，|ΔU_total/U_total| < 1e-10

7.2 守恒律断言的价值

守恒律是数学真理，任何正确实现都必然满足。如果实现中"通量双用"机制写错（例如 tank 用了被更新后的 flux，或算 flux 时没快照 W），守恒律立刻以可观察的数量级被破坏。这比"和参考解对比"更可靠——后者有过拟合风险。

7.3 不做的测试

• 绘图/动画输出（肉眼验证更高效）
• main.py 装配逻辑（无独立逻辑）
• 网格收敛 / CFL 收敛 / 阶数验证
• Method of Manufactured Solutions
• mock/fixture 框架（10 条测试不需要抽象）

7.4 运行方式

pytest tests/ -v

项目根运行。conftest.py 自动注入 src/ 到 sys.path。

7.5 可测性要求（对生产代码的约束）

1. 各模块顶层无副作用（禁止 print/read-file/random-seed）
2. solver.run() 不读取 config 模块，所有参数通过入参传入
3. Pipe.primitives() 返回 (ρ, u, P, a) 四元组，方便测试断言

---

8. 验收标准

MVP 完成的判据：

1. ✅ python src/main.py 能跑完、不崩溃、在 results/ 下生成全部 4 个文件
2. ✅ pytest tests/ -v 全部 10 条测试通过
3. ✅ main.py 末尾的"总质量 & 总能量相对误差 < 1e-10" sanity check 通过
4. ✅ P1(t) 单调递减、P2(t) 单调递增
5. ✅ 动画能用常见 GIF 播放器（浏览器、图片查看器）正常打开
6. ✅ history.npz 可以用 np.load 重新加载并访问所有字段
7. ✅ 全程运行时间（不含测试） < 30 秒

---

9. 已知局限与后续扩展

本 MVP 的物理局限：

• 细管（D=5mm）下摩擦损失可能显著，当前绝热无摩擦会高估压力传递速率
• 瓶内气体从滞止到出口的加速不满足严格等熵（因为 u_ghost=0 的简化）
• 一阶空间导致膨胀波被耗散
• 一阶时间精度对大 CFL 下的波传播有额外耗散

后续可扩展方向（不在本 MVP 范围内）：

1. 加 Darcy–Weisbach 摩擦源项（src/friction.py）
2. 加壁面对流换热（src/heat_transfer.py）
3. 升级到 MUSCL + minmod 二阶空间
4. 升级到 SSP-RK2 二阶时间
5. 在 Tank.ghost_state() 中加入入口损失系数 ζ
6. 配置化工况：YAML/JSON 读入 cases/*.yaml
7. 扩展到多管网络（带分叉节点的 0D/1D 混合拓扑）
8. 真实气体状态方程（Peng–Robinson / GERG-2008）

---

10. 参考

• Toro, E. F. Riemann Solvers and Numerical Methods for Fluid Dynamics. Springer. （HLL 格式的标准参考）
• Harten, A., Lax, P. D., and van Leer, B. (1983). "On upstream differencing and Godunov-type schemes for hyperbolic conservation laws." SIAM Review 25(1): 35–61.
• LeVeque, R. J. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press.
• 用户在 brainstorming 阶段提供的参考资料（涵盖 0D–1D 耦合的虚网格 + Riemann 方法）