# comac-cfd-coupling

**Repository Path**: wangyue2006/comac-cfd-coupling

## Basic Information

- **Project Name**: comac-cfd-coupling
- **Description**: 2026 上海开源大赛·AI+工业软件赛道·中国商飞命题【流体 CFD 方向】一三维耦合求解时间步长协调和边界耦合——方法验证原型(Python+numpy 自研一维喷嘴/轴对称 2.5D Euler 求解器与双向耦合框架, E1~E5 实验, 12 pytest)
- **Primary Language**: Unknown
- **License**: Apache-2.0
- **Default Branch**: main
- **Homepage**: None
- **GVP Project**: No

## Statistics

- **Stars**: 0
- **Forks**: 0
- **Created**: 2026-10-12
- **Last Updated**: 2026-10-12

## Categories & Tags

**Categories**: Uncategorized

**Tags**: None

## README

# comac-cfd-coupling

> 仓库地址：<https://gitee.com/wangyue2006/comac-cfd-coupling>

**2026 上海开源大赛 · AI+工业软件赛道 · 中国商飞命题**
**【流体 CFD 方向】一三维耦合求解时间步长协调和边界耦合 —— 方法验证原型**

[![tests](https://img.shields.io/badge/pytest-12%2F12-brightgreen)]() [![license](https://img.shields.io/badge/license-Apache--2.0-blue)]()

## 这是什么

针对商飞命题的三点技术要求，用 Python + numpy 自研**一维喷嘴求解器**、**轴对称(2.5D)简化三维求解器**与**一三维双向耦合框架**，在笔记本上完整复现并回答命题的核心问题：

1. **时间步长差异大** —— 一维侧 CFL 步长 ~10⁻⁹ s 量级、三维侧 ~10⁻⁶ s 量级（步长比 10:1 ~ 1000:1），用**子循环（subcycling）+ 子步采样线性插值同步**协调；
2. **边界耦合** —— 一维喷嘴出口 (流量 ṁ, 总温 T0, 静压 p, 密度 ρ) 以**通量注入式边界**传给三维进口；三维腔体**容积平均压力作为背压实时反馈**给一维出口边界；
3. **双向耦合精度与稳定性** —— 显式（零阶保持/线性外推）与隐式（固定欠松弛 / Aitken Δ² 迭代至界面残差收敛）两族耦合方案，在 10:1~1000:1 步长比下做了精度/稳定性/守恒性定量对比，并集成**毫秒级 PID 背压调节阀**演示控制逻辑与耦合步的同步。

**定位（诚实声明）**：这是一个**数值方法验证原型**，不是工业级 CFD。三维侧是轴对称(2.5D)降维模型、无粘 Euler 方程、无湍流/无相变；所有数字均由本仓库代码在 Mac 上真实运行产生。与工业级的差距见下文"原型边界"。

## 物理场景（缩比演示台）

空气（理想气体，p0=220 kPa, T0=300 K）经**收敛喷嘴**（L=10 mm，出口直径 40 mm）进入**等径稳压管段**（L=200 mm, R=20 mm，轴对称 2.5D 粗网格 20×3），出口为**简化阀门-孔口背压模型**通大气。初始全场静止 190 kPa（pr=0.86 > 临界 0.528，喷嘴全程亚声速，双向耦合始终有效），t=0 气源投入后充腔至平衡 ~206 kPa（管内 M≈0.35），并呈现无粘欠阻尼的腔-阀 Helmholtz 型振荡——为控制实验提供了真实动态。该拓扑与机舱供气/环控管路同构。

> 构型说明：曾尝试突扩腔（R=60mm）构型，其进口射流剪切层在无粘粗网格上存在离散 Kelvin–Helmholtz 不稳定（径向速度指数增长直至发散，诊断记录见实验报告附录）；工业求解器以湍流模型/人工耗散处理该问题。本原型改用等径构型，对耦合方法学（步长协调/边界耦合/稳定性）的验证无影响。

## 方法

### 一维求解器（`solver1d.py`）
准一维 Euler 方程有限体积：MUSCL(minmod) 重构 + Rusanov 通量 + Heun(RK2)；
面积加权守恒形式（均匀流严格保持）；亚声速进口给定 (p0, T0)、速度内场外推；
出口给定背压（壅塞时自动切换为超声速外推——此时背压物理上不起作用）。

### 三维求解器（`solver3d.py`）
轴对称 (x, r) Euler 方程有限体积（"2.5D"，柱坐标精确度量：面面积含环量因子，
含 r 动量几何源项 p·ΔA_r/V）：MUSCL + Rusanov + RK2；
进口面为"喷嘴口径开口 + 其余滑移壁"；出口为简化阀门模型
p_out = p_atm + ½kρu²/σ²（σ 为阀门开度）；步长采用**多维谱半径 CFL**
（含轴线单元 2/dr 几何刚性）。

### 耦合框架（`coupling.py`）—— 本题核心
```
正向:  1D 出口 (ṁ, T0, p, ρ) ──通量注入──> 3D 进口边界
反向:  3D 腔体容积平均静压 p_b ──背压边界──> 1D 出口
时间:  耦合步 dt_c = R·dt_1D, 窗口内一维子循环 R 步,
       出口参数按子步采样为时间序列, 三维步内按 RK2 阶段时刻线性插值取值
```
- **显式耦合** `explicit-hold` / `explicit-extrap`：背压用上一窗口值（零阶保持）或按趋势线性外推；一维子循环一次通过，无迭代；
- **隐式耦合** `implicit-fixed(ω)` / `implicit-aitken`：窗口内对耦合变量 p_b 迭代（Aitken Δ² 动态松弛 / 固定欠松弛），直至 |p_b^3D响应 − p_b^给定| < tol，对任意步长比无条件稳定，收敛后与"整体同步求解"一致；
- **界面守恒**：进口按"质量流量入口"推广——质量/动量/能量通量全部由一维出口状态确定（F = [q, q·uₑ+pₑ, 0, q·h₀]），界面输运量两侧严格一致，残余误差仅来自时间插值同步（实测 ~10⁻⁴ 量级）；
- **控制同步**（`controller.py`）：PID（抗积分饱和、测量微分）每 k 个耦合步采样一次腔压并输出阀门开度，与耦合数据交换严格同步。

## 快速开始

```bash
cd comac-cfd-coupling
python3 -m pytest tests/                 # 12 项测试（守恒/插值/极端步长比/PID）
PYTHONPATH=src python3 -m comac_coupling.experiments.run_e1_verification   # 求解器验证
PYTHONPATH=src python3 -m comac_coupling.experiments.run_e2_coupled_vs_uncoupled
PYTHONPATH=src python3 -m comac_coupling.experiments.run_e3_ratio_sweep     # 核心实验 ~30 min（含 R=1 参考解 30 万耦合窗口）
PYTHONPATH=src python3 -m comac_coupling.experiments.run_e4_relaxation
PYTHONPATH=src python3 -m comac_coupling.experiments.run_e5_control
# 或一次性: PYTHONPATH=src python3 -m comac_coupling.experiments.run_all
```
依赖：`numpy scipy matplotlib pytest`（Python ≥ 3.10）。结果输出至 `results/`（CSV/JSON/PNG），`实验报告.md` 中的全部数字即由这些脚本产生。

## 主要结果（真实运行数据，复现命令见上）

完整表格与图见 [实验报告.md](实验报告.md) 与 `results/`。摘要：

| 实验 | 结论 |
|---|---|
| E1 求解器验证 | 一维稳态喷嘴 vs 等熵解析：ṁ 误差 ≤1.8e-7（4 位以上有效数字一致），内部连续性 ρuA=const 偏差 ≤3.7e-5；光滑声脉冲收敛阶 ≈1.43（minmod TVD 典型行为）；轴对称与一维两个独立求解器交叉验证：光滑脉冲 L2(p)=1.9e-3、激波管 1.2e-2（二阶）；三维侧全局守恒达机器精度（<1e-10） |
| E2 耦合 vs 非耦合 | 非耦合（固定设计背压）的界面背压预测误差峰值 **25.5 kPa（运行点的 11%）**，流量时均偏差 **22.9%**；双向耦合界面质量/能量守恒误差 4e-5 |
| E3 步长比扫描（核心） | 10:1~1000:1 全部稳定（含显式）：L2(p_b) ≤1.8e-4；L2(ṁ) 随 R 增大（7e-5→9e-3，滞后+三维时间积分主导）；界面守恒误差 ∝ dt_c（1.1e-6→2.7e-4）；R=1000 时零阶保持比线性外推/隐式差 60-70%；三维侧**多维谱半径 CFL** 是 1000:1 的必要条件 |
| E4 松弛策略 | 界面扰动(+5 kPa)恢复：Aitken 3 次迭代（残差 4864→0.04 Pa）；ω=1.0 需 2 次；ω=0.1 迭代 25 次不收敛——小松弛因子并非更稳 |
| E5 控制同步 | PID 背压调节与耦合步同拍采样（采样周期 0.4~200 µs 扫描）：投入后 0.62 ms 压制腔-阀自然振荡并进入 2% 带，末态误差 ≤0.19%（快采样）/1.25%（200 µs 慢采样，整定时间 +45%） |

## 与命题三点要求的对照

| 命题要求 | 本仓库实现 | 状态 |
|---|---|---|
| 时间步长差异大（一维 µs 级 / 三维 ms 级）的协调 | 子循环 + 子步采样插值同步，步长比 10:1~1000:1 全覆盖（原型器件为厘米级缩比，绝对步长相应缩小，步长比量级与命题一致） | ✅ 已实现并定量验证 |
| 一维喷嘴瞬态流量/温度/相态参数向三维实时传递 | (ṁ, T0, p, ρ) 瞬态序列传递 + 通量注入边界；相态（单相气体假设）未实现 | ✅ 流量/温度已实现；⚠️ 相态为 TODO |
| 三维背压对一维的反向作用，保证双向耦合精度 | 腔压→一维背压边界，显式/隐式两族方案精度-稳定性定量对比，界面守恒 ~10⁻⁴ | ✅ 已实现 |
| 毫秒级控制逻辑同步 | PID 调节阀与耦合步同步采样，采样周期扫描 | ✅ 已实现 |

## 原型边界（与工业级的差距，诚实说明）

1. **三维侧是轴对称 2.5D**，非全 3D（无周向/ swirl/非轴对称几何）；无粘 Euler，无湍流模型、无壁面边界层、无传热；
2. **绝对时间尺度**为原型器件缩比（喷嘴 10 mm），一维步长 ~10⁻⁹ s、三维 ~10⁻⁶ s；工业场景（米级管网+腔体、三维隐式推进）对应 µs/ms 级，耦合算法（子循环/插值/迭代格式）完全相同；
3. **显式时间推进**（双侧 RK2+CFL）：工业三维 CFD 普遍隐式（SIMPLE/双时间步），步长不受声学 CFL 限制；本原型以网格尺度差（1000:1）复现等价的步长比；
4. **阀门-孔口为简化代数模型**（无堵塞修正、无阀门动态特性）；壅塞工况会物理性切断背压反馈（前向解耦），本原型选亚声速运行区间以保持双向耦合，壅塞切换逻辑在 `solver1d.py` 中实现但未作为主实验工况；
5. **相态传递未实现**（命题提到"相态参数"）：当前单组分理想气体；多相/相变界面传递为 TODO；
6. 未做网格无关性收敛研究（三维侧网格固定 30×9，一维侧 E1 做了收敛阶验证）。

## 目录结构

```
comac-cfd-coupling/
├── src/comac_coupling/
│   ├── gas.py           理想气体模型（两侧共用）
│   ├── solver1d.py      准一维喷嘴 Euler FV 求解器
│   ├── solver3d.py      轴对称(2.5D) Euler FV 求解器
│   ├── coupling.py      耦合驱动器（子循环/插值/显式/隐式Aitken）
│   ├── controller.py    PID 背压调节阀
│   ├── cases.py         算例定义
│   ├── metrics.py       误差/守恒指标
│   └── experiments/     E1-E5 实验脚本 + run_all
├── tests/               pytest 12 项
├── results/             实验输出（CSV/JSON/PNG，随仓库提交）
├── 实验报告.md          完整实验报告（框架/数字/结论/边界）
├── pyproject.toml  LICENSE (Apache-2.0)
```

## License

Apache-2.0（见 [LICENSE](LICENSE)）