# 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 数据结构

```python
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`

```python
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,
)
```

未来加载与访问示例：
```python
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 运行方式

```bash
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 方法）