merge: 并入水平卧式流化床空床流场仿真项目
# Conflicts: # .gitignore # README.md
This commit is contained in:
+2
-3
@@ -1,11 +1,10 @@
|
||||
# Python
|
||||
__pycache__/
|
||||
*.pyc
|
||||
*.pyo
|
||||
.venv/
|
||||
venv/
|
||||
|
||||
# 仿真输出（可由脚本重新生成）
|
||||
docs/assets/*.npz
|
||||
.ipynb_checkpoints/
|
||||
|
||||
# Node
|
||||
node_modules/
|
||||
|
||||
+118
@@ -0,0 +1,118 @@
|
||||
# PhysicsAiFluentAnalysis
|
||||
|
||||
> 物理世界流场实时仿真 —— 用 AI 神经算子把高保真 CFD 压缩到毫秒级推理。
|
||||
|
||||
[]()
|
||||
[]()
|
||||
[]()
|
||||
|
||||
---
|
||||
|
||||
## 一句话定义
|
||||
|
||||
以 **Ansys Fluent 高保真仿真为真值**，训练 **物理约束的 AI 代理模型（神经算子）**，在保持物理一致性的前提下，把小时级的流场计算压缩到 **毫秒级实时推理**，支撑数字孪生、工程设计优化与沉浸式可视化。
|
||||
|
||||
## 为什么这件事"炸裂"又"难"
|
||||
|
||||
| 维度 | 现状 | 目标 |
|
||||
| --- | --- | --- |
|
||||
| 单次 CFD 仿真 | 小时 ~ 天级 | 毫秒级（>1000x 加速） |
|
||||
| 参数扫描 | 一次一个工况 | 全参数空间连续可微 |
|
||||
| 交互性 | 离线出报告 | 实时拖拽参数、流场即时响应 |
|
||||
| 物理可信度 | 求解器保证 | 物理约束 + 守恒律 + 不确定性量化 |
|
||||
|
||||
**难在哪：** Navier-Stokes 方程的非线性与多尺度特性、湍流闭合问题、AI 长时推演的误差累积与发散、复杂几何的泛化、3D 大网格的实时推理算力瓶颈。这不是"套一个模型"能解决的，需要 **数据流水线 + 物理约束 + 工程加速 + 严格验证** 四线并进。
|
||||
|
||||
## 已落地 Demo：FFU 洁净室流场实时仿真
|
||||
|
||||
FFU（风机过滤单元）下方洁净区的流场、温湿度、压力与尘埃粒子输运 —— 同场景提供 **高保真真值** 与 **快速代理** 两条路径。
|
||||
|
||||
|
||||
|
||||
| | 高保真真值 | 快速代理 |
|
||||
| --- | --- | --- |
|
||||
| 网格 / 步数 | 360×270 / 14000 步 | 120×90 / 5000 步（同源求解器） |
|
||||
| 端到端耗时 | **76 s** | **7.3 s（16-27x 加速）** |
|
||||
| 趋势一致性 | 基准 | 速度 **r=0.95** · 粒子 r=0.83 · 温度 r=0.68 |
|
||||
|
||||
|
||||
|
||||
**主要趋势**：单向下扫排污染 → 绕障/缝隙加速（局部 ~1.5 m/s）→ 背风侧死角堆积滞留 → 送风正压/回风负压（< 2 Pa）。
|
||||
|
||||
详见 [docs/07-ffu-demo.md](docs/07-ffu-demo.md)。
|
||||
|
||||
## 系统总览
|
||||
|
||||
|
||||
|
||||
## 技术路线（一图流）
|
||||
|
||||
|
||||
|
||||
## 里程碑路线图
|
||||
|
||||
|
||||
|
||||
## 文档导航
|
||||
|
||||
| 文档 | 内容 |
|
||||
| --- | --- |
|
||||
| [docs/01-vision.md](docs/01-vision.md) | 愿景、痛点、目标与非目标、成功标准（KPI） |
|
||||
| [docs/02-architecture.md](docs/02-architecture.md) | 分层架构、数据流、技术选型与对比 |
|
||||
| [docs/03-roadmap.md](docs/03-roadmap.md) | 7 个阶段的目标、交付物与验收标准 |
|
||||
| [docs/04-data-and-ai.md](docs/04-data-and-ai.md) | Fluent 数据工厂、模型选型、训练与评估方案 |
|
||||
| [docs/05-risks.md](docs/05-risks.md) | 技术风险、缓解策略与开放问题 |
|
||||
| [docs/06-references.md](docs/06-references.md) | 论文、开源工具与数据集索引 |
|
||||
| [docs/07-ffu-demo.md](docs/07-ffu-demo.md) | FFU 洁净室 Demo：物理设置、趋势验证、复现命令 |
|
||||
|
||||
## 仓库结构
|
||||
|
||||
```
|
||||
PhysicsAiFluentAnalysis/
|
||||
├── README.md # 项目总览（本文件）
|
||||
├── package.json # 图件工具依赖 (SVG -> PNG)
|
||||
├── docs/ # 规划文档 + 演示图
|
||||
│ ├── 01-vision.md ~ 07-ffu-demo.md
|
||||
│ └── assets/ # PNG / SVG 演示图
|
||||
│ ├── architecture.* 系统架构图
|
||||
│ ├── pipeline.* 数据与训练流水线
|
||||
│ ├── roadmap.* 里程碑路线图
|
||||
│ ├── ffu-flow-field.* FFU 流场（高保真 LBM）
|
||||
│ ├── ffu-scalar-fields.png 温湿度/压力/粒子（高保真）
|
||||
│ ├── ffu-particle-tracks.png 粒子轨迹与滞留时间
|
||||
│ ├── ffu-flow-field-fast.png 快速代理流场（7 s）
|
||||
│ ├── ffu-scalar-fields-fast.png 快速代理标量场
|
||||
│ └── ffu-validation.png 趋势验证对比图
|
||||
└── tools/ # 可复现的生成脚本
|
||||
├── ffu_simulation.py # 高保真 D2Q9 LBM + 标量输运（76 s）
|
||||
├── ffu_surrogate.py # 快速代理（7 s）+ 趋势验证
|
||||
├── gen_diagrams.py # 架构/流水线/路线图 SVG 生成
|
||||
└── svg2png.js # SVG -> PNG 批量转换
|
||||
```
|
||||
|
||||
## 快速开始
|
||||
|
||||
```bash
|
||||
# 依赖: Python 3.11+ (numpy / numba / matplotlib), Node 18+ (仅图件转换需要)
|
||||
|
||||
# 1) 快速代理出图（约 7 秒，含趋势验证对比图）
|
||||
python tools/ffu_surrogate.py
|
||||
|
||||
# 2) 高保真真值（约 76 秒，输出 3 张图 + npz 场数据）
|
||||
python tools/ffu_simulation.py
|
||||
|
||||
# 3) 规划图（架构 / 流水线 / 路线图）
|
||||
python tools/gen_diagrams.py
|
||||
|
||||
# 4) SVG -> PNG
|
||||
npm install && node tools/svg2png.js
|
||||
```
|
||||
|
||||
## 当前状态
|
||||
|
||||
**P0 规划期**：愿景、架构、路线图、风险分析完成；FFU 洁净室 Demo 已可运行（高保真 + 快速代理双路径 + 趋势验证）。
|
||||
**下一步**：进入 P1，完成 2D 圆柱绕流 POC（Fluent 数据 → FNO 基线 → 误差评估）。
|
||||
|
||||
---
|
||||
|
||||
*本项目处于早期规划阶段，欢迎讨论与共建。*
|
||||
@@ -1,118 +1,222 @@
|
||||
# PhysicsAiFluentAnalysis
|
||||
# 水平卧式流化床 · 空床流场仿真
|
||||
|
||||
> 物理世界流场实时仿真 —— 用 AI 神经算子把高保真 CFD 压缩到毫秒级推理。
|
||||
|
||||
[]()
|
||||
[]()
|
||||
[]()
|
||||
> 用 GPU 加速的 **3D D3Q19 格子玻尔兹曼方法 (LBM)** 计算水平卧式流化床**空床**（无颗粒）气相流场，
|
||||
> 以**恒功率离心风机曲线**为压力进口，布风板为**多孔板**，流量与压差由物理自洽涌现。
|
||||
|
||||
---
|
||||
|
||||
## 一句话定义
|
||||
## 1. 物理场景
|
||||
|
||||
以 **Ansys Fluent 高保真仿真为真值**，训练 **物理约束的 AI 代理模型（神经算子）**，在保持物理一致性的前提下，把小时级的流场计算压缩到 **毫秒级实时推理**，支撑数字孪生、工程设计优化与沉浸式可视化。
|
||||
|
||||
## 为什么这件事"炸裂"又"难"
|
||||
|
||||
| 维度 | 现状 | 目标 |
|
||||
| --- | --- | --- |
|
||||
| 单次 CFD 仿真 | 小时 ~ 天级 | 毫秒级（>1000x 加速） |
|
||||
| 参数扫描 | 一次一个工况 | 全参数空间连续可微 |
|
||||
| 交互性 | 离线出报告 | 实时拖拽参数、流场即时响应 |
|
||||
| 物理可信度 | 求解器保证 | 物理约束 + 守恒律 + 不确定性量化 |
|
||||
|
||||
**难在哪：** Navier-Stokes 方程的非线性与多尺度特性、湍流闭合问题、AI 长时推演的误差累积与发散、复杂几何的泛化、3D 大网格的实时推理算力瓶颈。这不是"套一个模型"能解决的，需要 **数据流水线 + 物理约束 + 工程加速 + 严格验证** 四线并进。
|
||||
|
||||
## 已落地 Demo：FFU 洁净室流场实时仿真
|
||||
|
||||
FFU（风机过滤单元）下方洁净区的流场、温湿度、压力与尘埃粒子输运 —— 同场景提供 **高保真真值** 与 **快速代理** 两条路径。
|
||||
|
||||
|
||||
|
||||
| | 高保真真值 | 快速代理 |
|
||||
| --- | --- | --- |
|
||||
| 网格 / 步数 | 360×270 / 14000 步 | 120×90 / 5000 步（同源求解器） |
|
||||
| 端到端耗时 | **76 s** | **7.3 s（16-27x 加速）** |
|
||||
| 趋势一致性 | 基准 | 速度 **r=0.95** · 粒子 r=0.83 · 温度 r=0.68 |
|
||||
|
||||
|
||||
|
||||
**主要趋势**：单向下扫排污染 → 绕障/缝隙加速（局部 ~1.5 m/s）→ 背风侧死角堆积滞留 → 送风正压/回风负压（< 2 Pa）。
|
||||
|
||||
详见 [docs/07-ffu-demo.md](docs/07-ffu-demo.md)。
|
||||
|
||||
## 系统总览
|
||||
|
||||
|
||||
|
||||
## 技术路线（一图流）
|
||||
|
||||
|
||||
|
||||
## 里程碑路线图
|
||||
|
||||
|
||||
|
||||
## 文档导航
|
||||
|
||||
| 文档 | 内容 |
|
||||
| 项目 | 参数 |
|
||||
| --- | --- |
|
||||
| [docs/01-vision.md](docs/01-vision.md) | 愿景、痛点、目标与非目标、成功标准（KPI） |
|
||||
| [docs/02-architecture.md](docs/02-architecture.md) | 分层架构、数据流、技术选型与对比 |
|
||||
| [docs/03-roadmap.md](docs/03-roadmap.md) | 7 个阶段的目标、交付物与验收标准 |
|
||||
| [docs/04-data-and-ai.md](docs/04-data-and-ai.md) | Fluent 数据工厂、模型选型、训练与评估方案 |
|
||||
| [docs/05-risks.md](docs/05-risks.md) | 技术风险、缓解策略与开放问题 |
|
||||
| [docs/06-references.md](docs/06-references.md) | 论文、开源工具与数据集索引 |
|
||||
| [docs/07-ffu-demo.md](docs/07-ffu-demo.md) | FFU 洁净室 Demo：物理设置、趋势验证、复现命令 |
|
||||
| 布风板尺寸 | 2400 mm (x, 长度) × 567 mm (y, 宽度)，水平放置 |
|
||||
| 布风板厚度 | **2 mm** |
|
||||
| 布风板开孔 | 孔径 2 mm，开孔率 **5.6%**，孔轴 **5° 倾斜**（朝 -x，出料方向） |
|
||||
| 进风分区 | 沿长度均分 **4 段**（每段 600 mm），每段下方独立风箱 |
|
||||
| 风箱高度 | 450 mm |
|
||||
| 进风管 | 每段 **DN300**，从 567 mm 侧面水平居中进入（风向平行床板） |
|
||||
| 沸腾/渐扩段 | 由 2400 × 567 在 **1000 mm** 高度内渐扩至 2400 × 1200 |
|
||||
| 上部稳压段 | 2400 × 1200 × 400 mm |
|
||||
| 渐缩段 | 500 mm |
|
||||
| 排风 | 顶部居中 **DN800** 圆管 |
|
||||
|
||||
## 仓库结构
|
||||
### 工况（离心风机）
|
||||
|
||||
```
|
||||
PhysicsAiFluentAnalysis/
|
||||
├── README.md # 项目总览（本文件）
|
||||
├── package.json # 图件工具依赖 (SVG -> PNG)
|
||||
├── docs/ # 规划文档 + 演示图
|
||||
│ ├── 01-vision.md ~ 07-ffu-demo.md
|
||||
│ └── assets/ # PNG / SVG 演示图
|
||||
│ ├── architecture.* 系统架构图
|
||||
│ ├── pipeline.* 数据与训练流水线
|
||||
│ ├── roadmap.* 里程碑路线图
|
||||
│ ├── ffu-flow-field.* FFU 流场（高保真 LBM）
|
||||
│ ├── ffu-scalar-fields.png 温湿度/压力/粒子（高保真）
|
||||
│ ├── ffu-particle-tracks.png 粒子轨迹与滞留时间
|
||||
│ ├── ffu-flow-field-fast.png 快速代理流场（7 s）
|
||||
│ ├── ffu-scalar-fields-fast.png 快速代理标量场
|
||||
│ └── ffu-validation.png 趋势验证对比图
|
||||
└── tools/ # 可复现的生成脚本
|
||||
├── ffu_simulation.py # 高保真 D2Q9 LBM + 标量输运（76 s）
|
||||
├── ffu_surrogate.py # 快速代理（7 s）+ 趋势验证
|
||||
├── gen_diagrams.py # 架构/流水线/路线图 SVG 生成
|
||||
└── svg2png.js # SVG -> PNG 批量转换
|
||||
```
|
||||
| 项目 | 数值 |
|
||||
| --- | --- |
|
||||
| 单段额定风量 | 2500 m³/h（4 段合计 10000 m³/h） |
|
||||
| 额定压差 | 2900 Pa（布风板上下） |
|
||||
| 风机特性 | **恒功率输出** `Δp · Q = P`，P = 2900 × (10000/3600) ≈ **8056 W** |
|
||||
| 关死压力 | 4350 Pa（用于 Q→0 截断） |
|
||||
|
||||
## 快速开始
|
||||
> **关键物理认识**：2500 m³/h 与 2900 Pa 是风机**额定点**上的两个值，不是两个独立约束。
|
||||
> 离心风机给的是**风机曲线**，实际工况点是「风机曲线 ∩ 系统阻力曲线」的**唯一交点**。
|
||||
> 因此仿真中**只施加风机曲线压力**，流量由系统阻力自然涌现。
|
||||
|
||||
---
|
||||
|
||||
## 2. 数值方法
|
||||
|
||||
- **求解器**：D3Q19 格子玻尔兹曼，BGK 碰撞 + 推流 + 半步反弹（bounce-back）固体边界
|
||||
- **加速**：NVIDIA GPU（`numba.cuda`），实测约 **120 MLUPS**（百万格点更新/秒）
|
||||
- **网格**：**100 格/m**，240 × 130 × 265 = **8,268,000 单元**（本机 8 GB 显存上限）
|
||||
- **进口**：**压力进口**（固定密度 = 风机曲线压力），速度零梯度自由
|
||||
- **出口**：**压力出口**（密度 = 1，速度自由）
|
||||
- **布风板**：固体薄板 + 规则开孔（1 格孔 / 间距 4 格 ≈ 6.25% 开孔率）
|
||||
- **固体壁**：无滑移
|
||||
|
||||
### 缩放
|
||||
|
||||
| 量 | 说明 |
|
||||
| --- | --- |
|
||||
| `U_LB` | 0.015（进口参考速度对应的格子速度） |
|
||||
| `SCALE_U` | `U_in / U_LB` |
|
||||
| `SCALE_P` | `1.2 × SCALE_U²`（格子压力 → Pa） |
|
||||
|
||||
---
|
||||
|
||||
## 3. 结果（100 格/m，2 mm 板厚）
|
||||
|
||||
| 指标 | 数值 |
|
||||
| --- | --- |
|
||||
| 工况点风量 | **Q ≈ 1463 m³/h**（风机曲线与系统阻力交点） |
|
||||
| 风机压力 | 4350 Pa（关死） |
|
||||
| 风箱静压 | **4061 Pa** |
|
||||
| 沸腾段静压 | **181 Pa** |
|
||||
| **布风板压差** | **+3880 Pa**（风箱 > 沸腾段 ✓） |
|
||||
| 最大流速 | 4.89 m/s |
|
||||
| 求解 | 3900 步 / 184 s |
|
||||
|
||||
> 压力**符号、结构、量级**均正确：风箱高压 → 布风板节流降压 → 沸腾段低压 → 向排风递减。
|
||||
|
||||
---
|
||||
|
||||
## 4. 输出图（`docs/assets/`，19 张，4× 分辨率）
|
||||
|
||||
| 文件 | 内容 |
|
||||
| --- | --- |
|
||||
| `fbed-overview.png` | 四联总览（x-z / y-z / 近板 / 沸腾段） |
|
||||
| `fbed-montage-plenum.png` | 风箱层水平切面 ×5（4 段侧进风分布） |
|
||||
| `fbed-montage-nearplate.png` | 布风板近区水平切面 ×6（射流合并） |
|
||||
| `fbed-montage-freeboard.png` | 沸腾段/渐扩/排风水平切面 ×6 |
|
||||
| `fbed-montage-xz.png` | 垂直切面 x-z ×5 |
|
||||
| `fbed-montage-yz.png` | 垂直切面 y-z ×8（沿长度） |
|
||||
| `fbed-sections.png` | 4 个区段中心 y-z |
|
||||
| `fbed-zoom-sections.png` | 按区段剖分的 x-z 局部放大 |
|
||||
| `fbed-streamlines.png` | 双大图流线 |
|
||||
| `fbed-vectors.png` | 速度矢量图 |
|
||||
| `fbed-components.png` | 速度分量 ux/uy/uz |
|
||||
| `fbed-3d-cutaway.png` | 三正交切面三维剖切 |
|
||||
| `fbed-pressure.png` | 静压场四联 |
|
||||
| `fbed-pressure-drop.png` | 布风板压差分析 |
|
||||
| `fbed-uniformity.png` | 各区段流量均匀性 |
|
||||
| `fbed-decay.png` | 沿高度气速衰减 |
|
||||
| `fbed-flowrate.png` | 进风/布风板/排风 风量平衡 |
|
||||
| `fbed-tilt.png` | 布风板 5° 倾斜对近板风向影响 |
|
||||
| `fbed-geometry.png` | 几何与边界示意 |
|
||||
|
||||
数据：`docs/assets/fbed_fields.npz`（ux/uy/uz/rho/solid）。
|
||||
|
||||
---
|
||||
|
||||
## 5. 复现
|
||||
|
||||
```bash
|
||||
# 依赖: Python 3.11+ (numpy / numba / matplotlib), Node 18+ (仅图件转换需要)
|
||||
# 依赖: Python 3.11+, numpy, numba, numba-cuda, CUDA 13.2 runtime, matplotlib
|
||||
# 主仿真 (100 格/m, 约 3-5 分钟)
|
||||
python tools/fbed_simulation.py
|
||||
|
||||
# 1) 快速代理出图（约 7 秒，含趋势验证对比图）
|
||||
python tools/ffu_surrogate.py
|
||||
|
||||
# 2) 高保真真值（约 76 秒，输出 3 张图 + npz 场数据）
|
||||
python tools/ffu_simulation.py
|
||||
|
||||
# 3) 规划图（架构 / 流水线 / 路线图）
|
||||
python tools/gen_diagrams.py
|
||||
|
||||
# 4) SVG -> PNG
|
||||
npm install && node tools/svg2png.js
|
||||
# 批量出图 (秒级, 读 npz)
|
||||
python tools/fbed_figures.py
|
||||
python tools/fbed_analysis.py
|
||||
```
|
||||
|
||||
## 当前状态
|
||||
环境变量（可选）：
|
||||
|
||||
**P0 规划期**：愿景、架构、路线图、风险分析完成；FFU 洁净室 Demo 已可运行（高保真 + 快速代理双路径 + 趋势验证）。
|
||||
**下一步**：进入 P1，完成 2D 圆柱绕流 POC（Fluent 数据 → FNO 基线 → 误差评估）。
|
||||
| 变量 | 默认 | 说明 |
|
||||
| --- | --- | --- |
|
||||
| `FBED_PPM` | 100 | 网格分辨率 (格/m) |
|
||||
| `FBED_ULB` | 0.015 | 参考格子速度（越小越不可压，但越慢） |
|
||||
| `FBED_HOLE` / `FBED_SPACE` | 1 / 4 | 开孔大小 / 间距（格） |
|
||||
| `FBED_TB` | 120 | 求解墙钟预算 (s) |
|
||||
| `FBED_DP` | 2900 | 额定压差 (Pa) |
|
||||
|
||||
---
|
||||
|
||||
*本项目处于早期规划阶段，欢迎讨论与共建。*
|
||||
## 6. 脚本说明
|
||||
|
||||
| 文件 | 作用 |
|
||||
| --- | --- |
|
||||
| `tools/fbed_simulation.py` | 主求解器（GPU LBM + 风机曲线 + 多孔板） |
|
||||
| `tools/fbed_figures.py` | 切面/矢量/压力/流量等批量出图 |
|
||||
| `tools/fbed_analysis.py` | 几何/倾斜/均匀性/衰减/压差等分析图 |
|
||||
| `tools/fbed_check.py` / `fbed_check2.py` | 通量与压差诊断 |
|
||||
| `tools/fbed_probe.py` | 切面 ASCII 速度探针 |
|
||||
| `tools/fbed_rans.py` | RANS(投影法+AMG) 尝试版（未收敛，留档） |
|
||||
|
||||
---
|
||||
|
||||
## 7. 局限与已知问题（诚实声明）
|
||||
|
||||
1. **板厚亚网格**：2 mm 板厚在 100 格/m（格距 10 mm）下**小于一个网格**，只能按 1 格最薄板近似。要真正解析 2 mm 板 + 2 mm 孔，需 ~500 格/m（约 10.3 亿单元，需 ~157 GB 显存，本机 8 GB 无法承载）。
|
||||
2. **5° 倾角未施加**：倾角水平分量仅为垂直分量的 8.7%，粗网格下孔道倾斜不足 1 格，无法几何分辨；作为动量源施加时曾导致进口反流，暂关闭。方向已按设计设为 **-x**。
|
||||
3. **不可压近似**：LBM 密度波动 ~2.5%，压力场为相对值。
|
||||
4. **层流**：未加湍流模型，粗网格下湍流混合未解析。
|
||||
5. **工况点偏低**：细网格把布风板开孔阻力解析得更"阻"，工况点落在 1463 m³/h；若要标定到额定 10000 m³/h，需放大等效开孔或减薄板厚。
|
||||
|
||||
---
|
||||
|
||||
## 8. 前景价值与 AI 相关性
|
||||
|
||||
### 8.1 工业价值
|
||||
|
||||
流化床是**化工、能源、制药、食品**行业的核心装备（催化裂化、煤/生物质气化、干燥、造粒、燃烧、涂层）。
|
||||
**空床气相流场**是流化床设计与运行的基础：
|
||||
|
||||
- 布风均匀性 → 决定起始流化质量与死区分布
|
||||
- 孔板射流动量 → 决定气泡生成与颗粒夹带
|
||||
- 渐扩段气速衰减 → 决定颗粒扬析与返料
|
||||
- 风箱配气 → 决定分区控制精度
|
||||
|
||||
传统 CFD（Fluent / OpenFOAM）单次仿真为**小时～天级**，无法支撑设计迭代、参数扫描与在线优化。
|
||||
本项目用 **GPU LBM** 把空床流场压到**分钟级**，可直接用于：
|
||||
|
||||
| 场景 | 价值 |
|
||||
| --- | --- |
|
||||
| 布风板开孔/倾角优化 | 快速比选方案，缩短研发周期 |
|
||||
| 分区风量配比标定 | 验证各区段均匀性 |
|
||||
| 故障诊断（堵孔、偏流） | 快速复现异常工况 |
|
||||
| 教学/演示 | 实时可视化流场结构 |
|
||||
|
||||
### 8.2 科学意义
|
||||
|
||||
- 验证了 **「恒功率风机曲线 + 压力进口 + 多孔板」** 这一**物理自洽**的建模路线：
|
||||
流量与压差由「风机曲线 ∩ 系统阻力」唯一确定，而非人为指定。
|
||||
- 定量揭示了**不可压求解器表达高损失系统（损失系数 K≈900）的固有限制**：
|
||||
2900 Pa 压差在 LBM 中对应约 10% 密度波动，逼近不可压假设边界。
|
||||
- 给出了**亚网格多孔板的建模误差**（2 mm 板在 100 格/m 下不可分辨）。
|
||||
|
||||
### 8.3 AI 相关性（核心）
|
||||
|
||||
本项目天然是 **AI for Science / 神经算子** 的数据工厂与落地场景：
|
||||
|
||||
```
|
||||
GPU 求解器(真值) ──► 参数化数据集 ──► 神经算子训练 ──► 毫秒级实时推理
|
||||
(本仓库) (几何/边界/风速) (FNO/GNN/Transformer) (数字孪生)
|
||||
```
|
||||
|
||||
| AI 方向 | 与本项目的结合点 |
|
||||
| --- | --- |
|
||||
| **神经算子 (Neural Operator)** | FNO / DeepONet 把「解算子」参数化为网络权重；输入几何+边界，输出全场 |
|
||||
| **图神经网络 (GNN)** | MeshGraphNet 类模型适配非结构/复杂几何，天然处理不规则边界 |
|
||||
| **Transformer / Attention** | 把流场当序列建模，长程依赖捕捉回流与全局耦合 |
|
||||
| **物理约束机器学习 (PINN)** | 在损失中嵌入 Navier-Stokes / 守恒律，提升外推可信度 |
|
||||
| **数字孪生** | 毫秒推理 → 参数拖拽即时刷新 → 在线优化与预测性维护 |
|
||||
| **代理模型加速** | 用 1% 的算力逼近求解器精度，支撑大规模参数扫描 |
|
||||
|
||||
**闭环价值**：
|
||||
1. **数据**：GPU LBM 批量生成 (几何 + 风机曲线 + 开孔率) → 流场/压力场
|
||||
2. **训练**：FNO / GNN / Transformer 学习解算子
|
||||
3. **推理**：毫秒级出结果，替代逐次 CFD
|
||||
4. **迭代**：在线优化布风板与配气，反馈回设计
|
||||
|
||||
> 本仓库当前提供**第 1 环（真值引擎 + 快速基线）**；后续接入第 2–4 环即可实现
|
||||
> 「物理世界流场实时仿真」的完整闭环。
|
||||
|
||||
### 8.4 战略意义
|
||||
|
||||
- **自主可控**：全链路 Python + 开源栈（numba / CUDA / scipy / pyamg），无商业 CFD 许可依赖
|
||||
- **可扩展**：同一套 GPU 求解器可平移至洁净室、数据中心、燃烧室、气象等流场场景
|
||||
- **算力弹性**：单卡可跑 100 格/m，多卡/云 GPU 可扩展至 500 格/m 以上
|
||||
|
||||
---
|
||||
|
||||
## 9. 后续方向
|
||||
|
||||
- 用 **OpenFOAM**（`simpleFoam` + `k-ε` + `porousBafflePressure`）在 CPU 上跑更高分辨率（~250–300 格/m）
|
||||
- 布风板改用**多孔跳变边界**（Δp = K·½ρv²，K≈900），避免离散开孔的网格依赖
|
||||
- 以求解器生成参数化数据集 → 训练 **GNN/Transformer 神经算子** → 毫秒级实时推理
|
||||
|
||||
---
|
||||
|
||||
*本项目为物理世界流场实时仿真方向的水平卧式流化床 Demo。*
|
||||
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 1.0 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 2.7 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 415 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 406 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 554 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 2.0 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 1.8 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 1.5 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 1.8 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 1.8 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 3.4 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 899 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 623 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 1.7 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 2.9 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 1.6 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 716 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 7.2 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 1.3 MiB |
Binary file not shown.
@@ -0,0 +1,179 @@
|
||||
# 开发过程与提示词汇总
|
||||
|
||||
> 记录本项目从需求提出到最终交付的完整对话过程、关键技术决策与物理认识。
|
||||
> 项目：水平卧式流化床空床流场仿真（GPU LBM）
|
||||
|
||||
---
|
||||
|
||||
## 一、需求提出
|
||||
|
||||
**用户初始需求（原始提示词）：**
|
||||
|
||||
> 我的水平卧式流化床，床板开孔 2mm，开孔率 5.6%，床板的尺寸 2400x567mm，床板下方进风，
|
||||
> 进风分成 4 个区段，每一个区段都是一个箱体。每个区段都用直径 300 的进风管居中，从侧面进风，
|
||||
> 也就是进风风向和床板是平行的。床板是水平放，床板的孔是往出口方向倾斜了一点，不是正朝上的，
|
||||
> 倾斜很小，只有 5 度。进风风量是 2500 立方米每小时每个区段，然后进风风压，也就是床板上下方
|
||||
> 压差是 2900Pa，上面沸腾区段需要你模拟出来一个合适的放大形状，请模拟空床的流场，顶部沸腾段
|
||||
> 顶部是大排风管，还是上面这个快速计算的思路。
|
||||
|
||||
**需求拆解：**
|
||||
- 几何：2400×567 布风板、4 区段风箱、DN300 侧进风、2mm 孔 / 5.6% 开孔率 / 5° 倾斜、渐扩沸腾段、顶部大排风管
|
||||
- 工况：每段 2500 m³/h，压差 2900 Pa
|
||||
- 目标：空床气相流场，快速计算，结果可参考
|
||||
|
||||
---
|
||||
|
||||
## 二、关键澄清（问答）
|
||||
|
||||
用户对 3 个关键几何点做出选择：
|
||||
|
||||
| 问题 | 用户选择 |
|
||||
| --- | --- |
|
||||
| 4 段如何划分 / 进风管方向 | 沿长度均分 4×600mm，DN300 从 567 侧面水平居中进入 |
|
||||
| 沸腾段放大形状 | 矩形渐扩 + 顶部大排风（2.4×0.567 → 2.4×1.2，DN800） |
|
||||
| 孔板 5° 倾斜方向 | 朝 +x（x=2.4 为出料端）|
|
||||
|
||||
> **后记**：倾角方向在后续被用户更正为 **-x**（见"关键物理认识"第 4 条）。
|
||||
|
||||
---
|
||||
|
||||
## 三、技术路线演进
|
||||
|
||||
```
|
||||
LBM (GPU) ──► 发现问题 ──► RANS 尝试 ──► 失败 ──► 回到 LBM + 风机曲线
|
||||
```
|
||||
|
||||
### 阶段 1：GPU LBM 空床仿真
|
||||
- 3D D3Q19 LBM，`numba.cuda` 加速，实测 ~120 MLUPS
|
||||
- 布风板最初用「速度边界」强制表观气速
|
||||
- **问题**：压力场符号反了（顶部压力 > 风箱），且质量不守恒
|
||||
|
||||
### 阶段 2：RANS 尝试（用户选 A）
|
||||
- 用 `scipy` + `pyamg` 写投影法 + AMG 压力求解器
|
||||
- **失败**：共置网格棋盘振荡 + 进口速度边界过约束，压力爆到 -72703 Pa → NaN
|
||||
- 改交错网格 MAC，AMG 收敛正常（残差 9e-9），但速度进口仍把压力逼爆
|
||||
|
||||
### 阶段 3：回到 LBM + 用户的两个关键物理修正
|
||||
- **修正 1（恒功率风机曲线）**：`Δp · Q = P`，不再强制流量
|
||||
- **修正 2（压力进口）**：进口给风机压力，流量自然涌现
|
||||
- **结果**：压力符号彻底修正，风箱 > 沸腾段 ✓
|
||||
|
||||
---
|
||||
|
||||
## 四、关键物理认识（用户指出，扭转全局）
|
||||
|
||||
### 1. 不是"强制满足"，而是风机曲线
|
||||
> "2500 是风机额定流量，2900 是压力，并不是不可变的，你搞什么呢，不是强制满足，
|
||||
> 而是离心风机会出现一个风机曲线的"
|
||||
|
||||
**含义**：2500 m³/h 与 2900 Pa 是同一额定点上的两个值，**不是两个独立约束**。
|
||||
工况点是「风机曲线 ∩ 系统阻力曲线」的**唯一交点**。仿真中只施加风机曲线压力，流量涌现。
|
||||
|
||||
### 2. 恒功率输出
|
||||
> "风机是遵循恒功率的输出"
|
||||
|
||||
**含义**：风机曲线不是抛物线，而是 **`Δp · Q = P`（双曲线）**。
|
||||
额定点 (0.6944 m³/s, 2900 Pa) → P ≈ 2014 W/段，总 8056 W。
|
||||
|
||||
### 3. 不要纠结收敛
|
||||
> "你为什么总是在纠结于收敛，如果你走到收敛这条路，那么这个项目就不可能成功"
|
||||
|
||||
**含义**：目标是快速、物理自洽的结果，而非把残差压到机器精度。
|
||||
|
||||
### 4. 倾角方向
|
||||
> "孔的倾斜方向是这个绿色箭头方向，我觉得你搞错了"
|
||||
|
||||
**含义**：孔轴倾斜方向为 **上-左（-x）**，而非我之前设的 +x。
|
||||
|
||||
### 5. 板厚
|
||||
> "板厚只有 2mm，你按多少算的"
|
||||
|
||||
**含义**：布风板厚 **2 mm**，而我当时用了 50 mm（差 25 倍），已修正。
|
||||
|
||||
---
|
||||
|
||||
## 五、遇到的问题与修正
|
||||
|
||||
| # | 问题 | 根因 | 修正 |
|
||||
| --- | --- | --- | --- |
|
||||
| 1 | 压力场符号反了 | 布风板用速度边界，无压降；压力由几何决定 | 改为压力进口 + 多孔板 |
|
||||
| 2 | 质量不守恒（进 9903 / 出 20100） | 速度边界强制流量，与无散约束冲突 | 压力进口让流量自由 |
|
||||
| 3 | RANS 压力爆到 -72703 Pa | 共置网格棋盘 + 进口过约束 | 改交错网格；仍失败，放弃 |
|
||||
| 4 | 压力进口后流量反流 -5828 | 倾角动量源施加过强 | 暂关闭倾角力 |
|
||||
| 5 | 顶部"进风口"误判 | 渐扩段末端向排风管强汇流的贴壁高速带 | 探针证实为汇流，非源 |
|
||||
| 6 | 板厚 50mm | 我设错 | 改 2 mm |
|
||||
| 7 | 500 格/m 要求 | 硬件装不下 | 见"硬件限制" |
|
||||
|
||||
---
|
||||
|
||||
## 六、硬件限制（诚实核算）
|
||||
|
||||
| 网格 | 单元数 | D3Q19 双缓冲显存 |
|
||||
| --- | --- | --- |
|
||||
| 80 格/m | 4.2 M | 0.6 GB ✅ |
|
||||
| 100 格/m | 8.3 M | 1.3 GB ✅ |
|
||||
| 200 格/m | 66 M | 10.1 GB ❌ |
|
||||
| 300 格/m | 223 M | 33.9 GB ❌ |
|
||||
| 500 格/m | 1033 M | **157 GB** ❌ |
|
||||
|
||||
**本机**：RTX A2000 **8 GB 显存** + 64 GB 内存。
|
||||
- 500 格/m = 10.3 亿单元，光分布函数需 157 GB，**超显存 20 倍**，物理上装不下
|
||||
- 本机 GPU 实际上限 ≈ **100–150 格/m**
|
||||
- 要上 500 格/m，需云 GPU（A100/H100 80GB × 2–4）或 CPU 求解器（OpenFOAM，64GB 内存可跑 ~250–300 格/m）
|
||||
|
||||
---
|
||||
|
||||
## 七、最终交付
|
||||
|
||||
**仓库**：`http://124.220.63.9:30000/admin/Physics-Ai-Fluent-Analysis.git`（分支 `fluidized-bed`）
|
||||
|
||||
**内容**：
|
||||
- `tools/fbed_simulation.py`：主求解器（100 格/m，8.27M 单元）
|
||||
- `tools/fbed_figures.py` / `fbed_analysis.py`：出图
|
||||
- `docs/assets/`：19 张图（4× 分辨率）+ `fbed_fields.npz`
|
||||
- `README.md`：项目详细说明
|
||||
|
||||
**结果**：
|
||||
- 工况点 Q ≈ 1463 m³/h（风机曲线与系统阻力交点）
|
||||
- 风箱 4061 Pa > 沸腾段 181 Pa（**Δp = +3880 Pa，符号正确**）
|
||||
- 求解 3900 步 / 184 s
|
||||
|
||||
---
|
||||
|
||||
## 八、前景价值与 AI 相关性
|
||||
|
||||
### 工业价值
|
||||
流化床是化工/能源/制药/食品行业核心装备，空床流场决定布风均匀性、死区、颗粒夹带与配气精度。
|
||||
传统 CFD 小时~天级，无法支撑设计迭代；本 GPU LBM 把它压到**分钟级**，可用于布风板优化、分区配气标定、故障诊断与教学演示。
|
||||
|
||||
### 科学意义
|
||||
- 验证「恒功率风机曲线 + 压力进口 + 多孔板」物理自洽建模路线（工况点由两曲线交点唯一确定）
|
||||
- 定量揭示不可压求解器表达高损失系统（K≈900，2900 Pa）的固有限制
|
||||
- 给出亚网格多孔板的建模误差（2 mm 板在 100 格/m 下不可分辨）
|
||||
|
||||
### AI 相关性
|
||||
```
|
||||
GPU 求解器(真值) ──► 参数化数据集 ──► 神经算子训练 ──► 毫秒级实时推理
|
||||
(本仓库) (几何/边界/风速) (FNO/GNN/Transformer) (数字孪生)
|
||||
```
|
||||
| AI 方向 | 结合点 |
|
||||
| --- | --- |
|
||||
| 神经算子 FNO/DeepONet | 把解算子参数化为网络权重 |
|
||||
| GNN (MeshGraphNet) | 适配复杂/不规则几何 |
|
||||
| Transformer/Attention | 长程依赖，捕捉回流与全局耦合 |
|
||||
| PINN/物理约束 | 损失中嵌入 N-S 与守恒律 |
|
||||
| 数字孪生 | 毫秒推理 → 参数拖拽即时刷新 |
|
||||
|
||||
**战略意义**：全链路开源栈（numba/CUDA/scipy/pyamg），无商业 CFD 许可依赖；同一求解器可平移至洁净室、数据中心、燃烧室等流场场景。
|
||||
|
||||
---
|
||||
|
||||
## 九、后续方向
|
||||
|
||||
1. **OpenFOAM**：`simpleFoam` + `k-ε` + `porousBafflePressure`，CPU 跑 250–300 格/m
|
||||
2. **多孔跳变边界**替代离散开孔，消除网格依赖
|
||||
3. **神经算子**：以求解器生成参数化数据集 → 训练 GNN/Transformer → 毫秒级实时推理
|
||||
|
||||
---
|
||||
|
||||
*本文档汇总了从需求到交付的完整技术讨论，供后续迭代参考。*
|
||||
@@ -0,0 +1,400 @@
|
||||
#!/usr/bin/env python
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
水平卧式流化床 空床流场 —— 分析图 (几何/矢量/压差/流量/均匀性/衰减/倾斜/分段放大)
|
||||
读取 fbed_fields.npz。
|
||||
|
||||
python tools/fbed_analysis.py
|
||||
"""
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
import matplotlib
|
||||
matplotlib.use("Agg")
|
||||
import matplotlib.pyplot as plt
|
||||
from matplotlib.patches import Rectangle
|
||||
|
||||
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
||||
import fbed_simulation as fs
|
||||
|
||||
OUT = fs.OUT
|
||||
BG, FG, MUTED, ACCENT = fs.BG, fs.FG, fs.MUTED, fs.ACCENT
|
||||
CMAP_FLOW = fs.CMAP_FLOW
|
||||
XC, YC, ZC = fs.XC, fs.YC, fs.ZC
|
||||
X0, X1, Y0, Y1, Z0, Z1 = fs.X0, fs.X1, fs.Y0, fs.Y1, fs.Z0, fs.Z1
|
||||
DX, SCALE_U, BED_L, BED_W, N_SEC, SEC_L = fs.DX, fs.SCALE_U, fs.BED_L, fs.BED_W, fs.N_SEC, fs.SEC_L
|
||||
PPM = fs.PPM
|
||||
NZ, NY, NX = fs.NZ, fs.NY, fs.NX
|
||||
|
||||
D = np.load(os.path.join(OUT, "fbed_fields.npz"))
|
||||
UX, UY, UZ, RHO, SOLID = D["ux"], D["uy"], D["uz"], D["rho"], D["solid"]
|
||||
PLATE = D["plate"]
|
||||
FLUID = ~SOLID
|
||||
UX, UY, UZ = UX * SCALE_U, UY * SCALE_U, UZ * SCALE_U
|
||||
SPEED = np.sqrt(UX ** 2 + UY ** 2 + UZ ** 2)
|
||||
P = (RHO - 1.0) / 3.0 * (1.2 * SCALE_U ** 2)
|
||||
CELL_A = DX * DX
|
||||
IZP = fs._zi(0.0)
|
||||
|
||||
|
||||
def zi(v):
|
||||
return fs._zi(v)
|
||||
|
||||
|
||||
def yi(v):
|
||||
return fs._yi(v)
|
||||
|
||||
|
||||
def xi(v):
|
||||
return fs._xi(v)
|
||||
|
||||
|
||||
def quiver_slice(ax, which, idx, step=3, scale=None, title="", vmax=None):
|
||||
if which == "xz":
|
||||
ext = [X0, X1, Z0, Z1]; xl, yl = "x (m)", "z (m)"
|
||||
U, V, M = UX[:, idx, :], UZ[:, idx, :], SPEED[:, idx, :]
|
||||
m = SOLID[:, idx, :]
|
||||
xs, ys = XC, ZC
|
||||
elif which == "yz":
|
||||
ext = [Y0, Y1, Z0, Z1]; xl, yl = "y (m)", "z (m)"
|
||||
U, V, M = UY[:, :, idx], UZ[:, :, idx], SPEED[:, :, idx]
|
||||
m = SOLID[:, :, idx]
|
||||
xs, ys = YC, ZC
|
||||
else:
|
||||
ext = [X0, X1, Y0, Y1]; xl, yl = "x (m)", "y (m)"
|
||||
U, V, M = UX[idx, :, :], UY[idx, :, :], SPEED[idx, :, :]
|
||||
m = SOLID[idx, :, :]
|
||||
xs, ys = XC, YC
|
||||
U = np.where(m, np.nan, U); V = np.where(m, np.nan, V); M = np.where(m, np.nan, M)
|
||||
ax.imshow(np.where(m, np.nan, M), origin="lower", extent=ext, cmap=CMAP_FLOW,
|
||||
vmin=0, vmax=vmax or np.nanpercentile(M, 99), interpolation="bilinear",
|
||||
aspect="equal", alpha=0.85)
|
||||
Xs, Ys = np.meshgrid(xs[::step], ys[::step])
|
||||
Us = U[::step, ::step]; Vs = V[::step, ::step]; Ms = M[::step, ::step]
|
||||
if scale is None:
|
||||
scale = 6.0 * max(ext[1] - ext[0], ext[3] - ext[2]) / max(np.nanmax(np.hypot(Us, Vs)), 1e-6)
|
||||
ax.quiver(Xs, Ys, Us, Vs, Ms, cmap="cool", scale=scale, width=0.0022,
|
||||
clim=(0, np.nanpercentile(M, 99)))
|
||||
ax.contour(np.linspace(ext[0], ext[1], m.shape[1]),
|
||||
np.linspace(ext[2], ext[3], m.shape[0]),
|
||||
m.astype(float), levels=[0.5], colors="#adbac7", linewidths=0.8)
|
||||
ax.set_xlim(ext[0], ext[1]); ax.set_ylim(ext[2], ext[3])
|
||||
ax.set_xlabel(xl, fontsize=9); ax.set_ylabel(yl, fontsize=9)
|
||||
ax.tick_params(labelsize=8)
|
||||
if title:
|
||||
ax.set_title(title, fontsize=12)
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------
|
||||
def fig_geometry():
|
||||
fig, axes = plt.subplots(1, 2, figsize=(17.0, 7.6))
|
||||
for ax, which, idx in [(axes[0], "xz", yi(0.0)), (axes[1], "yz", xi(1.2))]:
|
||||
if which == "xz":
|
||||
ext = [X0, X1, Z0, Z1]; sol = SOLID[:, idx, :]
|
||||
xl, yl = "x (m)", "z (m)"
|
||||
else:
|
||||
ext = [Y0, Y1, Z0, Z1]; sol = SOLID[:, :, idx]
|
||||
xl, yl = "y (m)", "z (m)"
|
||||
ax.imshow(np.where(sol, 0.55, 0.08), origin="lower", extent=ext, cmap="gray",
|
||||
vmin=0, vmax=1, interpolation="nearest", aspect="equal")
|
||||
ax.contour(np.linspace(ext[0], ext[1], sol.shape[1]),
|
||||
np.linspace(ext[2], ext[3], sol.shape[0]),
|
||||
sol.astype(float), levels=[0.5], colors="#adbac7", linewidths=1.0)
|
||||
ax.set_xlabel(xl, fontsize=10); ax.set_ylabel(yl, fontsize=10)
|
||||
# 标注 (在 x-z 图)
|
||||
ax = axes[0]
|
||||
ax.annotate("4 x DN300 侧进风\n(平行床板, 2500 m³/h/段)", xy=(0.05, -0.225), xytext=(-0.35, 0.9),
|
||||
fontsize=10, color="#7ee787",
|
||||
arrowprops=dict(arrowstyle="->", color="#7ee787", lw=1.6))
|
||||
ax.annotate("布风板 2mm孔 5.6% 5°倾斜", xy=(1.2, 0.0), xytext=(1.55, -0.35),
|
||||
fontsize=10, color="#58a6ff",
|
||||
arrowprops=dict(arrowstyle="->", color="#58a6ff", lw=1.4))
|
||||
ax.annotate("顶部仅 DN800 排风\n(无进风)", xy=(1.2, 2.1), xytext=(1.45, 1.55),
|
||||
fontsize=10.5, color="#ff7b72",
|
||||
arrowprops=dict(arrowstyle="->", color="#ff7b72", lw=1.6))
|
||||
ax.set_title("几何 / 边界示意 (x-z, y=0)", fontsize=12)
|
||||
axes[1].set_title("几何 / 边界示意 (y-z, x=1.2)", fontsize=12)
|
||||
fig.suptitle("空床几何与边界条件：底部 4 段侧进风，顶部仅排风", fontsize=15, y=0.98)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.95])
|
||||
fig.savefig(os.path.join(OUT, "fbed-geometry.png"), dpi=540)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-geometry.png")
|
||||
|
||||
|
||||
def fig_vectors():
|
||||
fig, axes = plt.subplots(2, 2, figsize=(16.0, 12.5))
|
||||
quiver_slice(axes[0, 0], "xz", yi(0.0), step=3, title="(a) 垂直切面 y=0 速度矢量")
|
||||
quiver_slice(axes[0, 1], "yz", xi(1.2), step=3, title="(b) 垂直切面 x=1.2 速度矢量")
|
||||
quiver_slice(axes[1, 0], "xy", zi(0.15), step=3, title="(c) 布风板近区 z=0.15 速度矢量")
|
||||
quiver_slice(axes[1, 1], "xy", zi(0.90), step=3, title="(d) 沸腾段 z=0.90 速度矢量")
|
||||
fig.suptitle("空床流场速度矢量图（箭头方向=流向, 颜色=速度大小 m/s）", fontsize=15, y=0.99)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.97])
|
||||
fig.savefig(os.path.join(OUT, "fbed-vectors.png"), dpi=520)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-vectors.png")
|
||||
|
||||
|
||||
def fig_tilt():
|
||||
fig = plt.figure(figsize=(16.0, 10.5))
|
||||
gs = fig.add_gridspec(2, 2)
|
||||
ax = fig.add_subplot(gs[0, 0])
|
||||
iz = zi(0.05)
|
||||
quiver_slice(ax, "xy", iz, step=4, vmax=2.5,
|
||||
title="(a) 布风板正上方 z=0.05 矢量（放大看 +x 倾斜）")
|
||||
ax.set_xlim(0.6, 1.8); ax.set_ylim(-0.30, 0.30)
|
||||
|
||||
ax = fig.add_subplot(gs[0, 1])
|
||||
zs, ang, um = [], [], []
|
||||
for z in np.arange(0.03, 0.60, 0.02):
|
||||
iz = zi(z)
|
||||
m = FLUID[iz] & (np.abs(YC)[:, None] < BED_W / 2)
|
||||
if m.sum() < 10:
|
||||
continue
|
||||
ux_m, uz_m = UX[iz][m].mean(), UZ[iz][m].mean()
|
||||
zs.append(z); ang.append(np.degrees(np.arctan2(ux_m, uz_m))); um.append(ux_m)
|
||||
ax.plot(zs, ang, "-o", color=ACCENT, ms=4, label="实际偏角")
|
||||
ax.axhline(5.0, color="#ff7b72", ls="--", lw=1.2, label="理论 +5°")
|
||||
ax.axhline(0.0, color=MUTED, lw=0.8)
|
||||
ax.set_xlabel("高度 z (m)", fontsize=10); ax.set_ylabel("流向偏角 (°)", fontsize=10)
|
||||
ax.set_title("(b) 近板流向偏角随高度衰减", fontsize=12)
|
||||
ax.legend(fontsize=9); ax.grid(alpha=0.2)
|
||||
|
||||
ax = fig.add_subplot(gs[1, 0])
|
||||
zs2, uxs, uzs = [], [], []
|
||||
for z in np.arange(0.03, 0.60, 0.02):
|
||||
iz = zi(z)
|
||||
m = FLUID[iz] & (np.abs(YC)[:, None] < BED_W / 2)
|
||||
if m.sum() < 10:
|
||||
continue
|
||||
zs2.append(z); uxs.append(UX[iz][m].mean()); uzs.append(UZ[iz][m].mean())
|
||||
ax.plot(zs2, uzs, "-o", ms=4, color="#39d0d8", label="垂直分量 <uz>")
|
||||
ax.plot(zs2, uxs, "-s", ms=4, color="#ff7b72", label="水平分量 <ux>（5° 驱动）")
|
||||
ax.set_xlabel("高度 z (m)", fontsize=10); ax.set_ylabel("面平均速度 (m/s)", fontsize=10)
|
||||
ax.set_title("(c) 近板面平均速度分量", fontsize=12)
|
||||
ax.legend(fontsize=9); ax.grid(alpha=0.2)
|
||||
|
||||
ax = fig.add_subplot(gs[1, 1])
|
||||
xline = np.linspace(0.02, BED_L - 0.02, 60)
|
||||
iz = zi(0.06)
|
||||
prof = []
|
||||
for xv in xline:
|
||||
i = xi(xv)
|
||||
m = FLUID[iz, :, i]
|
||||
prof.append(UX[iz, :, i][m].mean() if m.sum() else np.nan)
|
||||
ax.plot(xline, prof, color="#ff7b72", lw=1.6)
|
||||
ax.axhline(0, color=MUTED, lw=0.8)
|
||||
ax.set_xlabel("x (m)", fontsize=10); ax.set_ylabel("<ux> (m/s)", fontsize=10)
|
||||
ax.set_title("(d) 近板水平速度沿长度分布（整体 +x 偏置）", fontsize=12)
|
||||
ax.grid(alpha=0.2)
|
||||
fig.suptitle("布风板 5° 倾斜对近板风向的影响", fontsize=15, y=0.99)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.97])
|
||||
fig.savefig(os.path.join(OUT, "fbed-tilt.png"), dpi=520)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-tilt.png")
|
||||
|
||||
|
||||
def fig_uniformity():
|
||||
fig, axes = plt.subplots(1, 3, figsize=(17.0, 5.4))
|
||||
qs = []
|
||||
for k in range(N_SEC):
|
||||
ix0, ix1 = xi(0.6 * k), xi(0.6 * (k + 1))
|
||||
q = UZ[IZP, :, ix0:ix1][FLUID[IZP, :, ix0:ix1]].sum() * CELL_A * 3600.0
|
||||
qs.append(q)
|
||||
ax = axes[0]
|
||||
bars = ax.bar([f"第{k+1}段" for k in range(N_SEC)], qs, color=ACCENT, alpha=0.85)
|
||||
ax.axhline(2500, color="#ff7b72", ls="--", lw=1.4, label="目标 2500 m³/h")
|
||||
for b, q in zip(bars, qs):
|
||||
ax.text(b.get_x() + b.get_width() / 2, q + 40, f"{q:.0f}", ha="center", fontsize=9)
|
||||
ax.set_ylabel("布风板风量 (m³/h)", fontsize=10); ax.legend(fontsize=9)
|
||||
ax.set_title("(a) 各区段风量", fontsize=12); ax.grid(alpha=0.2, axis="y")
|
||||
|
||||
ax = axes[1]
|
||||
prof = []
|
||||
xline = np.linspace(0.02, BED_L - 0.02, 120)
|
||||
for xv in xline:
|
||||
i = xi(xv)
|
||||
m = FLUID[IZP + 3, :, i]
|
||||
prof.append(UZ[IZP + 3, :, i][m].mean() if m.sum() else np.nan)
|
||||
ax.plot(xline, prof, color="#39d0d8", lw=1.8)
|
||||
ax.axhline(2500.0 / (BED_L * BED_W) / 3600.0 * 3600.0, color=MUTED, ls=":", lw=1)
|
||||
for k in range(1, N_SEC):
|
||||
ax.axvline(0.6 * k, color="#30363d", lw=0.8)
|
||||
ax.set_xlabel("x (m)", fontsize=10); ax.set_ylabel("近板 <uz> (m/s)", fontsize=10)
|
||||
ax.set_title("(b) 近板上升速度沿长度分布（均匀性）", fontsize=12); ax.grid(alpha=0.2)
|
||||
|
||||
ax = axes[2]
|
||||
yline = np.linspace(-BED_W / 2 + 0.02, BED_W / 2 - 0.02, 60)
|
||||
profy = []
|
||||
for yv in yline:
|
||||
j = yi(yv)
|
||||
m = FLUID[IZP + 3, j, :]
|
||||
profy.append(UZ[IZP + 3, j, :][m].mean() if m.sum() else np.nan)
|
||||
ax.plot(yline, profy, color="#7ee787", lw=1.8)
|
||||
ax.set_xlabel("y (m)", fontsize=10); ax.set_ylabel("近板 <uz> (m/s)", fontsize=10)
|
||||
ax.set_title("(c) 近板上升速度沿宽度分布（均匀性）", fontsize=12); ax.grid(alpha=0.2)
|
||||
fig.suptitle("布风板流量均匀性", fontsize=15, y=1.0)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.95])
|
||||
fig.savefig(os.path.join(OUT, "fbed-uniformity.png"), dpi=520)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-uniformity.png")
|
||||
|
||||
|
||||
def fig_decay():
|
||||
fig, axes = plt.subplots(1, 2, figsize=(15.0, 7.0))
|
||||
zs, qz, uz_avg, spd_avg, area = [], [], [], [], []
|
||||
for iz in range(NZ):
|
||||
m = FLUID[iz]
|
||||
if m.sum() < 20:
|
||||
continue
|
||||
zs.append(ZC[iz])
|
||||
qz.append(UZ[iz][m].sum() * CELL_A * 3600.0)
|
||||
uz_avg.append(UZ[iz][m].mean())
|
||||
spd_avg.append(SPEED[iz][m].mean())
|
||||
area.append(m.sum() * CELL_A)
|
||||
ax = axes[0]
|
||||
ax.plot(uz_avg, zs, "-o", ms=3, color="#39d0d8", label="面平均 <uz>")
|
||||
ax.plot(spd_avg, zs, "-s", ms=3, color="#ff7b72", label="面平均 |U|")
|
||||
ax.axvline(0, color=MUTED, lw=0.8)
|
||||
ax.axhline(0.0, color="#30363d", lw=0.8)
|
||||
ax.set_xlabel("速度 (m/s)", fontsize=10); ax.set_ylabel("高度 z (m)", fontsize=10)
|
||||
ax.set_title("(a) 沿高度气速衰减", fontsize=12); ax.legend(fontsize=9); ax.grid(alpha=0.2)
|
||||
|
||||
ax = axes[1]
|
||||
ax.plot(qz, zs, "-o", ms=3, color=ACCENT, label="通过水平面的净风量")
|
||||
ax.axvline(10000, color="#ff7b72", ls="--", lw=1.4, label="总风量 10000 m³/h")
|
||||
ax.set_xlabel("风量 (m³/h)", fontsize=10); ax.set_ylabel("高度 z (m)", fontsize=10)
|
||||
ax.set_title("(b) 沿高度风量守恒（连续性）", fontsize=12)
|
||||
ax.legend(fontsize=9); ax.grid(alpha=0.2)
|
||||
fig.suptitle("沿高度气速衰减与风量守恒", fontsize=15, y=0.99)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.96])
|
||||
fig.savefig(os.path.join(OUT, "fbed-decay.png"), dpi=520)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-decay.png")
|
||||
|
||||
|
||||
def fig_zoom_sections():
|
||||
fig, axes = plt.subplots(4, 1, figsize=(11.0, 16.0))
|
||||
for k, ax in enumerate(axes):
|
||||
quiver_slice(ax, "xz", yi(0.0), step=2, title=f"第 {k+1} 段 x∈[{0.6*k:.1f}, {0.6*(k+1):.1f}] m")
|
||||
ax.set_xlim(0.6 * k - 0.02, 0.6 * (k + 1) + 0.02)
|
||||
fig.suptitle("按区段剖分的 x-z 局部放大（含风箱进风/布风板/沸腾段）", fontsize=15, y=0.995)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.975])
|
||||
fig.savefig(os.path.join(OUT, "fbed-zoom-sections.png"), dpi=500)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-zoom-sections.png")
|
||||
|
||||
|
||||
def fig_pressure_drop():
|
||||
fig = plt.figure(figsize=(16.0, 10.5))
|
||||
gs = fig.add_gridspec(2, 2)
|
||||
ax = fig.add_subplot(gs[0, 0])
|
||||
ext = [X0, X1, Z0, Z1]
|
||||
Pxz = np.where(SOLID[:, yi(0), :], np.nan, P[:, yi(0), :])
|
||||
lim = np.nanpercentile(np.abs(Pxz), 99)
|
||||
im = ax.imshow(Pxz, origin="lower", extent=ext, cmap="RdBu_r", vmin=-lim, vmax=lim,
|
||||
aspect="equal", interpolation="bilinear")
|
||||
ax.axhline(0, color="#7ee787", lw=1.0, ls="--")
|
||||
ax.set_title("(a) 静压场 x-z (y=0)，绿线=布风板", fontsize=12)
|
||||
ax.set_xlabel("x (m)"); ax.set_ylabel("z (m)")
|
||||
cb = fig.colorbar(im, ax=ax, fraction=0.036, pad=0.02); cb.set_label("Pa", fontsize=9)
|
||||
|
||||
ax = fig.add_subplot(gs[0, 1])
|
||||
ext = [X0, X1, Y0, Y1]
|
||||
Pxy = np.where(SOLID[zi(0.10)], np.nan, P[zi(0.10)])
|
||||
im = ax.imshow(Pxy, origin="lower", extent=ext, cmap="RdBu_r",
|
||||
vmin=-np.nanpercentile(np.abs(Pxy), 99), vmax=np.nanpercentile(np.abs(Pxy), 99),
|
||||
aspect="equal", interpolation="bilinear")
|
||||
ax.set_title("(b) 布风板正上方 z=0.10 静压", fontsize=12)
|
||||
ax.set_xlabel("x (m)"); ax.set_ylabel("y (m)")
|
||||
cb = fig.colorbar(im, ax=ax, fraction=0.036, pad=0.02); cb.set_label("Pa", fontsize=9)
|
||||
|
||||
ax = fig.add_subplot(gs[1, 0])
|
||||
dps = []
|
||||
for k in range(N_SEC):
|
||||
ix0, ix1 = xi(0.6 * k), xi(0.6 * (k + 1))
|
||||
mp = FLUID[zi(-0.30):zi(-0.05), :, ix0:ix1]
|
||||
pp = P[zi(-0.30):zi(-0.05), :, ix0:ix1][mp].mean()
|
||||
mf = FLUID[zi(0.05):zi(0.20), :, ix0:ix1]
|
||||
pf = P[zi(0.05):zi(0.20), :, ix0:ix1][mf].mean()
|
||||
dps.append(pp - pf)
|
||||
bars = ax.bar([f"第{k+1}段" for k in range(N_SEC)], dps, color="#ff7b72", alpha=0.85)
|
||||
ax.axhline(2900, color="#7ee787", ls="--", lw=1.4, label="设计压差 2900 Pa")
|
||||
for b, v in zip(bars, dps):
|
||||
ax.text(b.get_x() + b.get_width() / 2, v + 5, f"{v:.0f}", ha="center", fontsize=9)
|
||||
ax.set_ylabel("布风板压差 (Pa)", fontsize=10); ax.legend(fontsize=9)
|
||||
ax.set_title("(c) 各段布风板压差（模型值）", fontsize=12); ax.grid(alpha=0.2, axis="y")
|
||||
|
||||
ax = fig.add_subplot(gs[1, 1])
|
||||
ext = [X0, X1, Z0, Z1]
|
||||
Pxz2 = np.where(SOLID[:, yi(0), :], np.nan, P[:, yi(0), :])
|
||||
ax.plot(XC, np.nanmean(P[:, yi(0), :], axis=0), color="#39d0d8", lw=1.4, label="沿 x 平均 (y=0)")
|
||||
for zc, col in [(zi(-0.30), "#58a6ff"), (zi(-0.05), "#7ee787"), (zi(0.10), "#ff7b72")]:
|
||||
ax.plot(XC, np.where(SOLID[zc, yi(0), :], np.nan, P[zc, yi(0), :]), lw=1.0,
|
||||
color=col, label=f"z={ZC[zc]:.2f} m")
|
||||
ax.set_xlabel("x (m)", fontsize=10); ax.set_ylabel("静压 (Pa)", fontsize=10)
|
||||
ax.set_title("(d) 静压沿长度剖面", fontsize=12); ax.legend(fontsize=8); ax.grid(alpha=0.2)
|
||||
fig.suptitle("压差分析：布风板上下压差与静压分布", fontsize=15, y=0.99)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.97])
|
||||
fig.savefig(os.path.join(OUT, "fbed-pressure-drop.png"), dpi=520)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-pressure-drop.png")
|
||||
|
||||
|
||||
def fig_flowrate():
|
||||
fig, axes = plt.subplots(1, 3, figsize=(17.0, 5.4))
|
||||
qs = []
|
||||
for k in range(N_SEC):
|
||||
ix0, ix1 = xi(0.6 * k), xi(0.6 * (k + 1))
|
||||
q = UZ[IZP, :, ix0:ix1][FLUID[IZP, :, ix0:ix1]].sum() * CELL_A * 3600.0
|
||||
qs.append(q)
|
||||
ax = axes[0]
|
||||
ax.bar([f"第{k+1}段" for k in range(N_SEC)], qs, color="#39d0d8", alpha=0.85)
|
||||
ax.axhline(2500, color="#ff7b72", ls="--", lw=1.4, label="目标 2500")
|
||||
ax.set_ylabel("风量 (m³/h)", fontsize=10); ax.legend(fontsize=9)
|
||||
ax.set_title("(a) 各区段布风板风量", fontsize=12); ax.grid(alpha=0.2, axis="y")
|
||||
|
||||
ax = axes[1]
|
||||
cum = np.cumsum([UZ[IZP, :, i][FLUID[IZP, :, i]].sum() for i in range(NX)]) * CELL_A * 3600.0
|
||||
ax.plot(XC, cum, color=ACCENT, lw=1.8)
|
||||
ax.axhline(10000, color="#ff7b72", ls="--", lw=1.2, label="总计 10000")
|
||||
ax.set_xlabel("x (m)", fontsize=10); ax.set_ylabel("累计风量 (m³/h)", fontsize=10)
|
||||
ax.set_title("(b) 沿长度累计风量", fontsize=12); ax.legend(fontsize=9); ax.grid(alpha=0.2)
|
||||
|
||||
ax = axes[2]
|
||||
iz_in = yi(fs.PIPE_Y0)
|
||||
qin = UY[:, iz_in, :][FLUID[:, iz_in, :]].sum() * CELL_A * 3600.0
|
||||
iz_out = NZ - 2
|
||||
qout = UZ[iz_out][FLUID[iz_out]].sum() * CELL_A * 3600.0
|
||||
qplate = sum(qs)
|
||||
labels = ["进风\n(DN300×4)", "布风板\n(通过)", "排风\n(DN800)"]
|
||||
vals = [qin, qplate, qout]
|
||||
cols = ["#7ee787", "#39d0d8", "#ff7b72"]
|
||||
b = ax.bar(labels, vals, color=cols, alpha=0.85)
|
||||
for bb, v in zip(b, vals):
|
||||
ax.text(bb.get_x() + bb.get_width() / 2, v + 60, f"{v:.0f}", ha="center", fontsize=9.5)
|
||||
ax.axhline(10000, color=MUTED, ls=":", lw=1)
|
||||
ax.set_ylabel("风量 (m³/h)", fontsize=10)
|
||||
ax.set_title("(c) 进风 / 布风板 / 排风 风量平衡", fontsize=12); ax.grid(alpha=0.2, axis="y")
|
||||
fig.suptitle("流量分析", fontsize=15, y=1.0)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.95])
|
||||
fig.savefig(os.path.join(OUT, "fbed-flowrate.png"), dpi=520)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-flowrate.png")
|
||||
|
||||
|
||||
def main():
|
||||
fig_geometry()
|
||||
fig_vectors()
|
||||
fig_tilt()
|
||||
fig_uniformity()
|
||||
fig_decay()
|
||||
fig_zoom_sections()
|
||||
fig_pressure_drop()
|
||||
fig_flowrate()
|
||||
print("[analysis] done")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -0,0 +1,58 @@
|
||||
#!/usr/bin/env python
|
||||
# -*- coding: utf-8 -*-
|
||||
"""空床流场诊断: 检查顶部是否有非物理进流 / 几何开口, 并采样 5 度倾斜效果。"""
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
||||
import fbed_simulation as fs
|
||||
|
||||
D = np.load(os.path.join(fs.OUT, "fbed_fields.npz"))
|
||||
UX, UY, UZ, SOLID = D["ux"], D["uy"], D["uz"], D["solid"]
|
||||
UX, UY, UZ = UX * fs.SCALE_U, UY * fs.SCALE_U, UZ * fs.SCALE_U
|
||||
fluid = ~SOLID
|
||||
NZ, NY, NX = SOLID.shape
|
||||
|
||||
print("=== 出口层 (z 顶部) 速度 ===")
|
||||
for iz in (NZ - 1, NZ - 2):
|
||||
m = fluid[iz]
|
||||
u = UZ[iz][m]
|
||||
print(f" iz={iz} z={fs.ZC[iz]:.3f} n={m.sum()} uz[min,mean,max]="
|
||||
f"[{u.min():.2f}, {u.mean():.2f}, {u.max():.2f}] m/s "
|
||||
f"inflow(<0)占比={100*np.mean(u<0):.1f}%")
|
||||
|
||||
print("\n=== 沸腾段顶部区域 (z>1.35) 是否有非物理入流 ===")
|
||||
band = fluid & (fs.ZC[:, None, None] > 1.35)
|
||||
uzb = UZ[band]
|
||||
print(f" n={uzb.size} uz mean={uzb.mean():.2f} inflow<0 占比={100*np.mean(uzb<0):.2f}%")
|
||||
|
||||
print("\n=== 采样疑似'进风口'位置 (yz 切面左上亮斑) ===")
|
||||
for (x, y, z) in [(0.15, -0.45, 1.2), (0.15, -0.35, 1.1), (2.25, -0.45, 1.2),
|
||||
(0.15, -0.50, 0.9), (1.05, -0.45, 1.2)]:
|
||||
ix, iy, iz = fs._xi(x), fs._yi(y), fs._zi(z)
|
||||
if fluid[iz, iy, ix]:
|
||||
print(f" (x={x},y={y},z={z}): |U|={np.hypot(np.hypot(UX[iz,iy,ix],UY[iz,iy,ix]),UZ[iz,iy,ix]):.2f} "
|
||||
f"ux={UX[iz,iy,ix]:+.2f} uy={UY[iz,iy,ix]:+.2f} uz={UZ[iz,iy,ix]:+.2f}")
|
||||
else:
|
||||
print(f" (x={x},y={y},z={z}): SOLID")
|
||||
|
||||
print("\n=== 5 度倾斜效果: 布风板正上方薄层水平速度 (y 平均) ===")
|
||||
for z in (0.04, 0.06, 0.10, 0.15):
|
||||
iz = fs._zi(z)
|
||||
m = fluid[iz] & (np.abs(fs.YC)[:, None] < fs.BED_W / 2)
|
||||
ux_m = UX[iz][m].mean()
|
||||
uz_m = UZ[iz][m].mean()
|
||||
print(f" z={z:.2f} <ux>={ux_m:+.3f} <uz>={uz_m:+.3f} m/s "
|
||||
f"偏角={np.degrees(np.arctan2(ux_m, uz_m)):+.2f} deg (理论 +5 deg)")
|
||||
|
||||
print("\n=== 布风板各段流量 (按 x 分 4 段) ===")
|
||||
izp = fs._zi(0.0)
|
||||
cell_a = fs.DX ** 2
|
||||
for k in range(4):
|
||||
x0, x1 = 0.6 * k, 0.6 * (k + 1)
|
||||
ix0, ix1 = fs._xi(x0), fs._xi(x1)
|
||||
sub = UZ[izp, :, ix0:ix1] * (fluid[izp, :, ix0:ix1])
|
||||
q = sub.sum() * cell_a * 3600.0
|
||||
print(f" 第{k+1}段 x[{x0:.1f},{x1:.1f}]: Q={q:.0f} m3/h")
|
||||
@@ -0,0 +1,50 @@
|
||||
#!/usr/bin/env python
|
||||
# -*- coding: utf-8 -*-
|
||||
"""沿高度净通量核查 (体积 / 质量)。"""
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
||||
import fbed_simulation as fs
|
||||
|
||||
D = np.load(os.path.join(fs.OUT, "fbed_fields.npz"))
|
||||
UX, UY, UZ, RHO, SOLID = D["ux"], D["uy"], D["uz"], D["rho"], D["solid"]
|
||||
FLUID = ~SOLID
|
||||
UX, UY, UZ = UX * fs.SCALE_U, UY * fs.SCALE_U, UZ * fs.SCALE_U
|
||||
CA = fs.DX ** 2
|
||||
|
||||
print(" z(m) 面积(m2) Q_vol(m3/h) Q_mass(m3/h) 倒流占比")
|
||||
for z in [-0.44, -0.35, -0.28, -0.22, -0.15, -0.08, -0.02, 0.02, 0.06, 0.12,
|
||||
0.25, 0.50, 0.80, 1.00, 1.20, 1.40, 1.60, 1.80, 1.95, 2.05, 2.15]:
|
||||
iz = fs._zi(z)
|
||||
m = FLUID[iz]
|
||||
if m.sum() < 5:
|
||||
continue
|
||||
area = m.sum() * CA
|
||||
qv = UZ[iz][m].sum() * CA * 3600
|
||||
qm = (RHO[iz] * UZ[iz])[m].sum() * CA * 3600
|
||||
frac = 100 * np.mean(UZ[iz][m] < 0)
|
||||
print(f" {z:+.2f} {area:7.3f} {qv:11.0f} {qm:11.0f} {frac:6.1f}%")
|
||||
|
||||
print("\n进/出总通量:")
|
||||
iy_in = fs._yi(fs.PIPE_Y0)
|
||||
m = FLUID[:, iy_in, :]
|
||||
print(f" 进风口(4xDN300) Q_vol={UY[:, iy_in, :][m].sum()*CA*3600:.0f} "
|
||||
f"Q_mass={(RHO*UY)[:, iy_in, :][m].sum()*CA*3600:.0f} m3/h")
|
||||
for iz in (fs.NZ - 1, fs.NZ - 2):
|
||||
m = FLUID[iz]
|
||||
print(f" 排风口 z={fs.ZC[iz]:.3f} Q_vol={UZ[iz][m].sum()*CA*3600:.0f} "
|
||||
f"Q_mass={(RHO[iz]*UZ[iz])[m].sum()*CA*3600:.0f} m3/h")
|
||||
|
||||
print("\n=== 压差方向核查 (布风板上下) ===")
|
||||
P = (RHO - 1.0) / 3.0 * (1.2 * fs.SCALE_U ** 2)
|
||||
for k in range(4):
|
||||
ix0, ix1 = fs._xi(0.6 * k), fs._xi(0.6 * (k + 1))
|
||||
mpl = FLUID[fs._zi(-0.35):fs._zi(-0.05), :, ix0:ix1]
|
||||
p_pl = P[fs._zi(-0.35):fs._zi(-0.05), :, ix0:ix1][mpl].mean()
|
||||
mfb = FLUID[fs._zi(0.05):fs._zi(0.20), :, ix0:ix1]
|
||||
p_fb = P[fs._zi(0.05):fs._zi(0.20), :, ix0:ix1][mfb].mean()
|
||||
print(f" 第{k+1}段: p_风箱={p_pl:7.1f} Pa p_沸腾段={p_fb:7.1f} Pa "
|
||||
f"Δp(风箱-沸腾段)={p_pl-p_fb:+7.1f} Pa {'✓正确' if p_pl>p_fb else '✗反了'}")
|
||||
@@ -0,0 +1,312 @@
|
||||
#!/usr/bin/env python
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
水平卧式流化床 空床流场 —— 批量出图
|
||||
====================================
|
||||
读取 fbed_fields.npz, 生成大量切面 / 分量 / 压力 / 分区图。
|
||||
|
||||
python tools/fbed_figures.py
|
||||
"""
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
import matplotlib
|
||||
matplotlib.use("Agg")
|
||||
import matplotlib.pyplot as plt
|
||||
from matplotlib.colors import LinearSegmentedColormap
|
||||
|
||||
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
||||
import fbed_simulation as fs
|
||||
|
||||
OUT = fs.OUT
|
||||
BG, FG, MUTED, ACCENT = fs.BG, fs.FG, fs.MUTED, fs.ACCENT
|
||||
CMAP_FLOW = fs.CMAP_FLOW
|
||||
XC, YC, ZC = fs.XC, fs.YC, fs.ZC
|
||||
X0, X1, Y0, Y1, Z0, Z1 = fs.X0, fs.X1, fs.Y0, fs.Y1, fs.Z0, fs.Z1
|
||||
PPM, DX, SCALE_U = fs.PPM, fs.DX, fs.SCALE_U
|
||||
|
||||
D = np.load(os.path.join(OUT, "fbed_fields.npz"))
|
||||
UX, UY, UZ, RHO, SOLID = D["ux"], D["uy"], D["uz"], D["rho"], D["solid"]
|
||||
SPEED = np.sqrt(UX ** 2 + UY ** 2 + UZ ** 2) * SCALE_U
|
||||
P = (RHO - 1.0) / 3.0 * (1.2 * SCALE_U ** 2)
|
||||
P[SOLID] = np.nan
|
||||
|
||||
CMAP_P = plt.get_cmap("RdBu_r")
|
||||
|
||||
|
||||
def xi(v):
|
||||
return fs._xi(v)
|
||||
|
||||
|
||||
def yi(v):
|
||||
return fs._yi(v)
|
||||
|
||||
|
||||
def zi(v):
|
||||
return fs._zi(v)
|
||||
|
||||
|
||||
def get(which, idx):
|
||||
if which == "xz":
|
||||
return (SPEED[:, idx, :], UX[:, idx, :], UZ[:, idx, :], SOLID[:, idx, :],
|
||||
[X0, X1, Z0, Z1], "x (m)", "z (m)")
|
||||
if which == "yz":
|
||||
return (SPEED[:, :, idx], UY[:, :, idx], UZ[:, :, idx], SOLID[:, :, idx],
|
||||
[Y0, Y1, Z0, Z1], "y (m)", "z (m)")
|
||||
return (SPEED[idx, :, :], UX[idx, :, :], UY[idx, :, :], SOLID[idx, :, :],
|
||||
[X0, X1, Y0, Y1], "x (m)", "y (m)")
|
||||
|
||||
|
||||
def get_scalar(field, which, idx):
|
||||
if which == "xz":
|
||||
return field[:, idx, :], SOLID[:, idx, :], [X0, X1, Z0, Z1], "x (m)", "z (m)"
|
||||
if which == "yz":
|
||||
return field[:, :, idx], SOLID[:, :, idx], [Y0, Y1, Z0, Z1], "y (m)", "z (m)"
|
||||
return field[idx, :, :], SOLID[idx, :, :], [X0, X1, Y0, Y1], "x (m)", "y (m)"
|
||||
|
||||
|
||||
def panel(ax, which, idx, vmax, stream=True, density=1.1):
|
||||
s2, u2, v2, sol2, ext, xl, yl = get(which, idx)
|
||||
ax.imshow(np.where(sol2, np.nan, s2), origin="lower", extent=ext, cmap=CMAP_FLOW,
|
||||
vmin=0, vmax=vmax, interpolation="bilinear", aspect="equal")
|
||||
if stream:
|
||||
ny, nx = sol2.shape
|
||||
xs = np.linspace(ext[0], ext[1], nx)
|
||||
ys = np.linspace(ext[2], ext[3], ny)
|
||||
U = np.where(sol2, np.nan, u2)
|
||||
V = np.where(sol2, np.nan, v2)
|
||||
with np.errstate(invalid="ignore"):
|
||||
ax.streamplot(xs, ys, U, V, density=density,
|
||||
color=(0.85, 0.88, 0.92, 0.38), linewidth=0.45,
|
||||
arrowsize=0.5, zorder=5)
|
||||
ax.contour(np.linspace(ext[0], ext[1], sol2.shape[1]),
|
||||
np.linspace(ext[2], ext[3], sol2.shape[0]),
|
||||
sol2.astype(float), levels=[0.5], colors="#adbac7", linewidths=0.7, zorder=6)
|
||||
ax.set_xlabel(xl, fontsize=8)
|
||||
ax.set_ylabel(yl, fontsize=8)
|
||||
ax.set_xlim(ext[0], ext[1])
|
||||
ax.set_ylim(ext[2], ext[3])
|
||||
ax.tick_params(labelsize=7)
|
||||
|
||||
|
||||
def panel_scalar(ax, field, which, idx, cmap, vmin=None, vmax=None, center=False, stream=True):
|
||||
f2, sol2, ext, xl, yl = get_scalar(field, which, idx)
|
||||
f2 = np.where(sol2, np.nan, f2)
|
||||
if center:
|
||||
lim = np.nanpercentile(np.abs(f2), 99) * 1.05
|
||||
im = ax.imshow(f2, origin="lower", extent=ext, cmap=cmap, vmin=-lim, vmax=lim,
|
||||
interpolation="bilinear", aspect="equal")
|
||||
else:
|
||||
im = ax.imshow(f2, origin="lower", extent=ext, cmap=cmap, vmin=vmin, vmax=vmax,
|
||||
interpolation="bilinear", aspect="equal")
|
||||
if stream:
|
||||
s2, u2, v2, _, _, _, _ = get(which, idx)
|
||||
ny, nx = sol2.shape
|
||||
xs = np.linspace(ext[0], ext[1], nx)
|
||||
ys = np.linspace(ext[2], ext[3], ny)
|
||||
with np.errstate(invalid="ignore"):
|
||||
ax.streamplot(xs, ys, np.where(sol2, np.nan, u2), np.where(sol2, np.nan, v2),
|
||||
density=1.0, color=(0.85, 0.88, 0.92, 0.3), linewidth=0.4,
|
||||
arrowsize=0.5, zorder=5)
|
||||
ax.contour(np.linspace(ext[0], ext[1], sol2.shape[1]),
|
||||
np.linspace(ext[2], ext[3], sol2.shape[0]),
|
||||
sol2.astype(float), levels=[0.5], colors="#adbac7", linewidths=0.7, zorder=6)
|
||||
ax.set_xlabel(xl, fontsize=8)
|
||||
ax.set_ylabel(yl, fontsize=8)
|
||||
ax.set_xlim(ext[0], ext[1])
|
||||
ax.set_ylim(ext[2], ext[3])
|
||||
ax.tick_params(labelsize=7)
|
||||
return im
|
||||
|
||||
|
||||
def montage(specs, ncols, suptitle, name, vmax, cbar_label="|U| (m/s)",
|
||||
figw=17.0, panel_h=3.0, stream=True, density=1.1):
|
||||
nrows = (len(specs) + ncols - 1) // ncols
|
||||
fig, axes = plt.subplots(nrows, ncols, figsize=(figw, panel_h * nrows + 1.0),
|
||||
squeeze=False)
|
||||
for ax in axes.ravel()[len(specs):]:
|
||||
ax.axis("off")
|
||||
for ax, (spec, title) in zip(axes.ravel(), specs):
|
||||
which, idx = spec
|
||||
panel(ax, which, idx, vmax, stream=stream, density=density)
|
||||
ax.set_title(title, fontsize=9.5, pad=4)
|
||||
sm = plt.cm.ScalarMappable(cmap=CMAP_FLOW, norm=plt.Normalize(0, vmax))
|
||||
cb = fig.colorbar(sm, ax=axes.ravel().tolist(), fraction=0.012, pad=0.01)
|
||||
cb.set_label(cbar_label, fontsize=10)
|
||||
cb.outline.set_edgecolor("#30363d")
|
||||
fig.suptitle(suptitle, fontsize=15, y=0.995)
|
||||
fig.savefig(os.path.join(OUT, name), dpi=500, bbox_inches="tight")
|
||||
plt.close(fig)
|
||||
print(" ->", name)
|
||||
|
||||
|
||||
def fig_xy_plenum():
|
||||
specs = [(("xy", zi(z)), f"z = {z:+.2f} m") for z in (-0.42, -0.32, -0.22, -0.12, -0.02)]
|
||||
montage(specs, 3, "风箱层水平切面 |U|（4 段侧进风分布）", "fbed-montage-plenum.png", 6.0)
|
||||
|
||||
|
||||
def fig_xy_plate():
|
||||
specs = [(("xy", zi(z)), f"z = {z:.2f} m") for z in (0.04, 0.08, 0.13, 0.20, 0.30, 0.45)]
|
||||
montage(specs, 3, "布风板近区水平切面 |U|（射流合并过程）", "fbed-montage-nearplate.png", 3.0)
|
||||
|
||||
|
||||
def fig_xy_free():
|
||||
specs = [(("xy", zi(z)), f"z = {z:.2f} m") for z in (0.60, 0.80, 1.00, 1.20, 1.45, 1.75)]
|
||||
montage(specs, 3, "沸腾段 / 渐扩 / 排风 水平切面 |U|", "fbed-montage-freeboard.png", 4.0)
|
||||
|
||||
|
||||
def fig_xz():
|
||||
specs = [(("xz", yi(y)), f"y = {y:+.2f} m") for y in (-0.26, -0.13, 0.0, 0.13, 0.26)]
|
||||
montage(specs, 3, "垂直切面 x-z（不同宽度位置）|U|", "fbed-montage-xz.png", 5.0, figw=16.0)
|
||||
|
||||
|
||||
def fig_yz():
|
||||
specs = [(("yz", xi(x)), f"x = {x:.2f} m") for x in (0.15, 0.45, 0.75, 1.05, 1.35, 1.65, 1.95, 2.25)]
|
||||
montage(specs, 4, "垂直切面 y-z（沿长度）|U|", "fbed-montage-yz.png", 5.0, figw=18.0, panel_h=3.4)
|
||||
|
||||
|
||||
def fig_components():
|
||||
fig, axes = plt.subplots(2, 2, figsize=(15.0, 12.0))
|
||||
panel(axes[0, 0], "xz", yi(0.0), 4.0)
|
||||
axes[0, 0].set_title("(a) 速度大小 |U|", fontsize=12, color=ACCENT)
|
||||
for ax, comp, title, cmap in [
|
||||
(axes[0, 1], UX, "(b) 水平分量 ux（孔板 5° 倾斜 + 向排风汇流）", "coolwarm"),
|
||||
(axes[1, 0], UY, "(c) 横向分量 uy（侧进风 + 渐扩）", "coolwarm"),
|
||||
(axes[1, 1], UZ, "(d) 垂直分量 uz（上升主流）", "coolwarm")]:
|
||||
f2, sol2, ext, xl, yl = get_scalar(comp * SCALE_U, "xz", yi(0.0))
|
||||
lim = np.nanpercentile(np.abs(np.where(sol2, np.nan, f2)), 99) * 1.05
|
||||
im = ax.imshow(np.where(sol2, np.nan, f2), origin="lower", extent=ext, cmap=cmap,
|
||||
vmin=-lim, vmax=lim, interpolation="bilinear", aspect="equal")
|
||||
ax.contour(np.linspace(ext[0], ext[1], sol2.shape[1]),
|
||||
np.linspace(ext[2], ext[3], sol2.shape[0]),
|
||||
sol2.astype(float), levels=[0.5], colors="#adbac7", linewidths=0.7)
|
||||
ax.set_title(title, fontsize=12)
|
||||
ax.set_xlabel("x (m)", fontsize=9)
|
||||
ax.set_ylabel("z (m)", fontsize=9)
|
||||
cb = fig.colorbar(im, ax=ax, fraction=0.036, pad=0.02)
|
||||
cb.set_label("m/s", fontsize=9)
|
||||
cb.outline.set_edgecolor("#30363d")
|
||||
fig.suptitle("速度分量分解（垂直切面 y=0）", fontsize=15, y=0.99)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.97])
|
||||
fig.savefig(os.path.join(OUT, "fbed-components.png"), dpi=520)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-components.png")
|
||||
|
||||
|
||||
def fig_pressure():
|
||||
fig, axes = plt.subplots(2, 2, figsize=(15.0, 12.0))
|
||||
for ax, (which, idx, title) in zip(axes.ravel(), [
|
||||
("xz", yi(0.0), "(a) 垂直切面 y=0 —— 全系统压差"),
|
||||
("yz", xi(1.2), "(b) 垂直切面 x=1.2"),
|
||||
("xy", zi(0.15), "(c) 布风板近区 z=0.15"),
|
||||
("xy", zi(0.90), "(d) 沸腾段 z=0.90")]):
|
||||
im = panel_scalar(ax, P, which, idx, CMAP_P, center=True, stream=False)
|
||||
ax.set_title(title, fontsize=12)
|
||||
cb = fig.colorbar(im, ax=ax, fraction=0.036, pad=0.02)
|
||||
cb.set_label("静压 (Pa)", fontsize=9)
|
||||
cb.outline.set_edgecolor("#30363d")
|
||||
fig.suptitle("静压场（相对出口, 单位 Pa）", fontsize=15, y=0.99)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.97])
|
||||
fig.savefig(os.path.join(OUT, "fbed-pressure.png"), dpi=520)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-pressure.png")
|
||||
|
||||
|
||||
def fig_sections():
|
||||
fig, axes = plt.subplots(1, 4, figsize=(17.0, 6.0))
|
||||
for k, ax in enumerate(axes):
|
||||
xc = 0.6 * (k + 0.5)
|
||||
panel(ax, "yz", xi(xc), 5.0)
|
||||
ax.set_title(f"第 {k+1} 段 x={xc:.1f} m", fontsize=11)
|
||||
fig.suptitle("4 个进风区段中心垂直切面 y-z（单段风量 2500 m³/h）", fontsize=14, y=1.0)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.95])
|
||||
fig.savefig(os.path.join(OUT, "fbed-sections.png"), dpi=520)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-sections.png")
|
||||
|
||||
|
||||
def fig_hero_stream():
|
||||
fig, axes = plt.subplots(1, 2, figsize=(18.0, 8.0))
|
||||
panel(axes[0], "xz", yi(0.0), 4.0, density=1.8)
|
||||
axes[0].set_title("垂直切面 y=0 —— 风箱/布风板/沸腾段/排风", fontsize=12)
|
||||
panel(axes[1], "yz", xi(1.2), 4.0, density=1.8)
|
||||
axes[1].set_title("垂直切面 x=1.2 —— 侧向渐扩", fontsize=12)
|
||||
fig.suptitle("空床流场流线（GPU D3Q19 LBM）", fontsize=15, y=0.99)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.96])
|
||||
fig.savefig(os.path.join(OUT, "fbed-streamlines.png"), dpi=560)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-streamlines.png")
|
||||
|
||||
|
||||
def fig_3d():
|
||||
from mpl_toolkits.mplot3d import Axes3D # noqa
|
||||
fig = plt.figure(figsize=(13.0, 10.0))
|
||||
ax = fig.add_subplot(111, projection="3d")
|
||||
norm = plt.Normalize(0, 4.0)
|
||||
|
||||
def surf_y(iy):
|
||||
f = SPEED[:, iy, :]
|
||||
m = SOLID[:, iy, :]
|
||||
X, Z = np.meshgrid(XC, ZC)
|
||||
Y = np.full_like(X, YC[iy])
|
||||
fc = CMAP_FLOW(norm(np.where(m, np.nan, f)))
|
||||
fc[m, 3] = 0.0
|
||||
ax.plot_surface(X, Y, Z, facecolors=fc, shade=False, rstride=1, cstride=1,
|
||||
linewidth=0, antialiased=False)
|
||||
|
||||
def surf_x(ix):
|
||||
f = SPEED[:, :, ix]
|
||||
m = SOLID[:, :, ix]
|
||||
Y, Z = np.meshgrid(YC, ZC)
|
||||
X = np.full_like(Y, XC[ix])
|
||||
fc = CMAP_FLOW(norm(np.where(m, np.nan, f)))
|
||||
fc[m, 3] = 0.0
|
||||
ax.plot_surface(X, Y, Z, facecolors=fc, shade=False, rstride=1, cstride=1,
|
||||
linewidth=0, antialiased=False)
|
||||
|
||||
def surf_z(iz):
|
||||
f = SPEED[iz, :, :]
|
||||
m = SOLID[iz, :, :]
|
||||
X, Y = np.meshgrid(XC, YC)
|
||||
Z = np.full_like(X, ZC[iz])
|
||||
fc = CMAP_FLOW(norm(np.where(m, np.nan, f)))
|
||||
fc[m, 3] = 0.0
|
||||
ax.plot_surface(X, Y, Z, facecolors=fc, shade=False, rstride=1, cstride=1,
|
||||
linewidth=0, antialiased=False)
|
||||
|
||||
surf_y(yi(0.0))
|
||||
surf_x(xi(1.2))
|
||||
surf_z(zi(0.9))
|
||||
ax.set_xlabel("x (m)")
|
||||
ax.set_ylabel("y (m)")
|
||||
ax.set_zlabel("z (m)")
|
||||
ax.set_title("三正交切面三维剖切（|U|）", fontsize=13)
|
||||
ax.view_init(elev=22, azim=-58)
|
||||
sm = plt.cm.ScalarMappable(cmap=CMAP_FLOW, norm=norm)
|
||||
cb = fig.colorbar(sm, ax=ax, fraction=0.03, pad=0.08, shrink=0.6)
|
||||
cb.set_label("|U| (m/s)", fontsize=10)
|
||||
cb.outline.set_edgecolor("#30363d")
|
||||
fig.savefig(os.path.join(OUT, "fbed-3d-cutaway.png"), dpi=520)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-3d-cutaway.png")
|
||||
|
||||
|
||||
def main():
|
||||
print("[figs] loading", os.path.join(OUT, "fbed_fields.npz"))
|
||||
fig_xy_plenum()
|
||||
fig_xy_plate()
|
||||
fig_xy_free()
|
||||
fig_xz()
|
||||
fig_yz()
|
||||
fig_components()
|
||||
fig_pressure()
|
||||
fig_sections()
|
||||
fig_hero_stream()
|
||||
fig_3d()
|
||||
print("[figs] done")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -0,0 +1,54 @@
|
||||
#!/usr/bin/env python
|
||||
# -*- coding: utf-8 -*-
|
||||
"""探针: 检查 x=0.3 的 y-z 切面是否有异常开口 / 高速射流。"""
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
||||
import fbed_simulation as fs
|
||||
|
||||
D = np.load(os.path.join(fs.OUT, "fbed_fields.npz"))
|
||||
UX, UY, UZ, SOLID = D["ux"], D["uy"], D["uz"], D["solid"]
|
||||
UX, UY, UZ = UX * fs.SCALE_U, UY * fs.SCALE_U, UZ * fs.SCALE_U
|
||||
SP = np.sqrt(UX ** 2 + UY ** 2 + UZ ** 2)
|
||||
NZ, NY, NX = SOLID.shape
|
||||
|
||||
ix = fs._xi(0.3)
|
||||
print(f"=== x=0.3 (ix={ix}) y-z 切面: 固体/流体 ASCII 图 (z 由高到低) ===")
|
||||
print("每格: '#'=固体 '.'=流体(<0.5m/s) 'o'=0.5-1.5 'O'=1.5-3 '*'=3-6 '@'=>6 m/s")
|
||||
for iz in range(NZ - 1, -1, -1):
|
||||
z = fs.ZC[iz]
|
||||
if z < -0.5 or z > 1.7:
|
||||
continue
|
||||
row = ""
|
||||
for iy in range(NY):
|
||||
if SOLID[iz, iy, ix]:
|
||||
row += "#"
|
||||
else:
|
||||
s = SP[iz, iy, ix]
|
||||
row += "." if s < 0.5 else ("o" if s < 1.5 else ("O" if s < 3 else ("*" if s < 6 else "@")))
|
||||
print(f"z={z:+.2f} {row}")
|
||||
|
||||
print("\n=== y=-0.6 墙完整性 (沿 z 扫) ===")
|
||||
iy_wall = fs._yi(-0.6)
|
||||
runs = []
|
||||
for iz in range(NZ):
|
||||
if SOLID[iz, iy_wall, ix]:
|
||||
runs.append((fs.ZC[iz], "solid"))
|
||||
else:
|
||||
runs.append((fs.ZC[iz], "FLUID"))
|
||||
open_z = [z for z, s in runs if s == "FLUID"]
|
||||
print(f" y=-0.6 处流体 z 段: {min(open_z):.2f} .. {max(open_z):.2f}" if open_z else " 全固体")
|
||||
print(f" 流体 z 计数={len(open_z)} / {NZ}")
|
||||
|
||||
print("\n=== 顶部区域 (z>1.0) 是否存在指向域内的射流 ===")
|
||||
for (x, y, z) in [(0.3, -0.5, 1.2), (0.3, -0.45, 1.3), (0.3, -0.35, 1.1),
|
||||
(2.1, -0.5, 1.2), (0.9, -0.5, 1.2), (1.5, -0.5, 1.2)]:
|
||||
i = fs._xi(x); j = fs._yi(y); k = fs._zi(z)
|
||||
if SOLID[k, j, i]:
|
||||
print(f" ({x},{y},{z}): SOLID")
|
||||
else:
|
||||
print(f" ({x},{y},{z}): ux={UX[k,j,i]:+6.2f} uy={UY[k,j,i]:+6.2f} uz={UZ[k,j,i]:+6.2f} "
|
||||
f"|U|={SP[k,j,i]:5.2f}")
|
||||
@@ -0,0 +1,401 @@
|
||||
#!/usr/bin/env python
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
水平卧式流化床 空床流场 —— RANS + 多孔跳变布风板 (交错网格投影法 + AMG)
|
||||
======================================================================
|
||||
交错网格 (MAC): p 在单元中心, u/v/w 在面心 -> 从原理上消除棋盘压力振荡。
|
||||
|
||||
预测: u* = u + dt[ -(u.grad)u + div(nu_eff grad u) + S_porous ]
|
||||
压力: div(grad p) = (rho/dt) div(u*) (pyamg 多重网格, 固体 Neumann, 出口 Dirichlet)
|
||||
修正: u = u* - (dt/rho) grad(p)
|
||||
|
||||
布风板 = 薄多孔层 (Darcy), 阻力标定 dP=2900 Pa @ 表观气速 2.04 m/s。
|
||||
进风 = 4 x DN300 速度进口 (2500 m3/h/段); 排风 = 压力出口; 其余无滑移壁。
|
||||
|
||||
运行:
|
||||
python tools/fbed_rans.py
|
||||
"""
|
||||
import os
|
||||
import time
|
||||
|
||||
import numpy as np
|
||||
import numba
|
||||
from scipy.sparse import csr_matrix
|
||||
import pyamg
|
||||
|
||||
# --------------------------------------------------------------------------
|
||||
BED_L, BED_W = 2.4, 0.567
|
||||
N_SEC = 4
|
||||
SEC_L = BED_L / N_SEC
|
||||
PLENUM_H = 0.45
|
||||
PLATE_T = 0.05
|
||||
FB_H, FB_W = 1.0, 1.2
|
||||
UPPER_H = 0.4
|
||||
REDUCER_H = 0.5
|
||||
DUCT_H = 0.3
|
||||
EXH_R = 0.4
|
||||
PIPE_R = 0.15
|
||||
PIPE_Y0 = -0.65
|
||||
|
||||
X0, X1 = 0.0, BED_L
|
||||
Y0, Y1 = PIPE_Y0, 0.65
|
||||
Z0 = -PLENUM_H
|
||||
Z1 = FB_H + UPPER_H + REDUCER_H + DUCT_H
|
||||
|
||||
PPM = int(os.environ.get("RANS_PPM", "20"))
|
||||
NX = int(round((X1 - X0) * PPM))
|
||||
NY = int(round((Y1 - Y0) * PPM))
|
||||
NZ = int(round((Z1 - Z0) * PPM))
|
||||
DX = (X1 - X0) / NX
|
||||
DY = (Y1 - Y0) / NY
|
||||
DZ = (Z1 - Z0) / NZ
|
||||
XC = X0 + (np.arange(NX) + 0.5) * DX
|
||||
YC = Y0 + (np.arange(NY) + 0.5) * DY
|
||||
ZC = Z0 + (np.arange(NZ) + 0.5) * DZ
|
||||
|
||||
Q_SEC = 2500.0 / 3600.0
|
||||
A_PIPE = np.pi * PIPE_R ** 2
|
||||
U_IN = Q_SEC / A_PIPE
|
||||
A_BED = BED_L * BED_W
|
||||
U_SUP = Q_SEC * N_SEC / A_BED
|
||||
DP_RATED = float(os.environ.get("RANS_DP", "2900.0")) # 额定点压力 (2500 m3/h/段)
|
||||
Q_RATED = Q_SEC * N_SEC # 额定总流量 m3/s
|
||||
P_FAN = DP_RATED * Q_RATED # 恒功率 (W)
|
||||
DP_SHUT = DP_RATED * float(os.environ.get("RANS_FAN_SHUT", "1.5"))
|
||||
|
||||
|
||||
def fan_dp(q_m3s):
|
||||
"""恒功率风机曲线: dp * Q = P -> dp = P/Q (Q->0 时用关死压力截断)"""
|
||||
if q_m3s < 1e-6:
|
||||
return DP_SHUT
|
||||
return min(P_FAN / q_m3s, DP_SHUT)
|
||||
RHO = 1.2
|
||||
NU = 1.5e-5
|
||||
R_POROUS = DP_RATED / (U_SUP * PLATE_T)
|
||||
NU_T = float(os.environ.get("RANS_NUT", "0.02"))
|
||||
NU_EFF = NU + NU_T
|
||||
|
||||
MAX_STEPS = int(os.environ.get("RANS_STEPS", "4000"))
|
||||
TIME_BUDGET = float(os.environ.get("RANS_TB", "16.0"))
|
||||
CFL = float(os.environ.get("RANS_CFL", "0.3"))
|
||||
|
||||
|
||||
def build_masks():
|
||||
solid = np.ones((NZ, NY, NX), dtype=bool)
|
||||
yw = BED_W / 2.0
|
||||
|
||||
def box(xa, xb, ya, yb, za, zb, val=False):
|
||||
ix0, ix1 = int(round((xa - X0) / DX)), int(round((xb - X0) / DX))
|
||||
iy0, iy1 = int(round((ya - Y0) / DY)), int(round((yb - Y0) / DY))
|
||||
iz0, iz1 = int(round((za - Z0) / DZ)), int(round((zb - Z0) / DZ))
|
||||
solid[max(iz0, 0):min(iz1, NZ), max(iy0, 0):min(iy1, NY),
|
||||
max(ix0, 0):min(ix1, NX)] = val
|
||||
|
||||
box(0.0, BED_L, -yw, yw, -PLENUM_H, 0.0)
|
||||
for k in range(1, N_SEC):
|
||||
xk = k * SEC_L
|
||||
box(xk - 0.0125, xk + 0.0125, -yw, yw, -PLENUM_H, 0.0, True)
|
||||
box(0.0, BED_L, -yw, yw, 0.0, PLATE_T)
|
||||
for iz in range(NZ):
|
||||
z = ZC[iz]
|
||||
if z <= PLATE_T or z >= FB_H:
|
||||
continue
|
||||
w = yw + (FB_W / 2.0 - yw) * (z / FB_H)
|
||||
solid[iz] = np.where(np.abs(YC[:, None]) <= w, False, solid[iz])
|
||||
box(0.0, BED_L, -FB_W / 2.0, FB_W / 2.0, FB_H, FB_H + UPPER_H)
|
||||
z_r0, z_r1 = FB_H + UPPER_H, FB_H + UPPER_H + REDUCER_H
|
||||
for iz in range(NZ):
|
||||
z = ZC[iz]
|
||||
if z < z_r0 or z >= z_r1:
|
||||
continue
|
||||
t = (z - z_r0) / REDUCER_H
|
||||
xh = BED_L / 2.0 - (BED_L / 2.0 - 0.4) * t
|
||||
yh = FB_W / 2.0 - (FB_W / 2.0 - 0.4) * t
|
||||
solid[iz] = np.where((np.abs(XC[None, :] - BED_L / 2.0) <= xh) &
|
||||
(np.abs(YC[:, None]) <= yh), False, solid[iz])
|
||||
cyl = (XC[None, :] - BED_L / 2.0) ** 2 + (YC[:, None]) ** 2 <= EXH_R ** 2
|
||||
for iz in range(NZ):
|
||||
if ZC[iz] >= z_r1:
|
||||
solid[iz] = np.where(cyl, False, solid[iz])
|
||||
|
||||
iy_p0 = int(round((PIPE_Y0 - Y0) / DY))
|
||||
iy_p1 = int(round((-yw - Y0) / DY))
|
||||
for k in range(N_SEC):
|
||||
xk = SEC_L * (k + 0.5)
|
||||
zk = -PLENUM_H / 2.0
|
||||
for iz in range(NZ):
|
||||
if abs(ZC[iz] - zk) > PIPE_R:
|
||||
continue
|
||||
m = (XC - xk) ** 2 + (ZC[iz] - zk) ** 2 <= PIPE_R ** 2
|
||||
solid[iz, iy_p0:iy_p1, m] = False
|
||||
|
||||
plate = np.zeros_like(solid)
|
||||
iz0, iz1 = int(round((0.0 - Z0) / DZ)), int(round((PLATE_T - Z0) / DZ))
|
||||
iy0, iy1 = int(round((-yw - Y0) / DY)), int(round((yw - Y0) / DY))
|
||||
plate[iz0:iz1, iy0:iy1, :] = ~solid[iz0:iz1, iy0:iy1, :]
|
||||
inlet = np.zeros_like(solid)
|
||||
inlet[:, iy_p0:iy_p0 + 2, :] = ~solid[:, iy_p0:iy_p0 + 2, :]
|
||||
outlet = np.zeros_like(solid)
|
||||
outlet[NZ - 2:, :, :] = ~solid[NZ - 2:, :, :]
|
||||
return solid, plate, inlet, outlet
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------
|
||||
# 预测步 (交错网格)
|
||||
# --------------------------------------------------------------------------
|
||||
@numba.njit(cache=True, fastmath=True)
|
||||
def predict(u, v, w, solid, plate, inlet, outlet, nu, dt, dx, dy, dz, rporous, uin):
|
||||
nz, ny, nx = solid.shape
|
||||
# u faces: (nz, ny, nx+1), u[iz,iy,ix] at west face of cell ix
|
||||
for iz in range(nz):
|
||||
for iy in range(ny):
|
||||
for ix in range(1, nx):
|
||||
if solid[iz, iy, ix] and solid[iz, iy, ix - 1]:
|
||||
u[iz, iy, ix] = 0.0
|
||||
continue
|
||||
if solid[iz, iy, ix] or solid[iz, iy, ix - 1]:
|
||||
u[iz, iy, ix] = 0.0
|
||||
continue
|
||||
uc = u[iz, iy, ix]
|
||||
ue = u[iz, iy, ix + 1] if ix + 1 <= nx else uc
|
||||
uw = u[iz, iy, ix - 1]
|
||||
un = u[iz, iy + 1, ix] if iy + 1 < ny else uc
|
||||
us = u[iz, iy - 1, ix] if iy - 1 >= 0 else uc
|
||||
ut = u[iz + 1, iy, ix] if iz + 1 < nz else uc
|
||||
ub = u[iz - 1, iy, ix] if iz - 1 >= 0 else uc
|
||||
vn = 0.5 * (v[iz, iy + 1, ix - 1] + v[iz, iy + 1, ix]) if iy + 1 < ny else 0.0
|
||||
vs = 0.5 * (v[iz, iy, ix - 1] + v[iz, iy, ix])
|
||||
wt = 0.5 * (w[iz + 1, iy, ix - 1] + w[iz + 1, iy, ix]) if iz + 1 < nz else 0.0
|
||||
wb = 0.5 * (w[iz, iy, ix - 1] + w[iz, iy, ix])
|
||||
dux = (uc - uw) / dx if uc > 0 else (ue - uc) / dx
|
||||
duy = (un - uc) / dy if vn > 0 else (uc - us) / dy
|
||||
duz = (ut - uc) / dz if wt > 0 else (uc - ub) / dz
|
||||
conv = uc * dux + 0.5 * (vn + vs) * duy + 0.5 * (wt + wb) * duz
|
||||
lap = (ue - 2.0 * uc + uw) / (dx * dx) + (un - 2.0 * uc + us) / (dy * dy) \
|
||||
+ (ut - 2.0 * uc + ub) / (dz * dz)
|
||||
unew = uc + dt * (-conv + nu * lap)
|
||||
if plate[iz, iy, ix] or plate[iz, iy, ix - 1]:
|
||||
unew /= (1.0 + dt * rporous)
|
||||
u[iz, iy, ix] = unew
|
||||
# v faces: (nz, ny+1, nx)
|
||||
for iz in range(nz):
|
||||
for iy in range(1, ny):
|
||||
for ix in range(nx):
|
||||
if solid[iz, iy, ix] or solid[iz, iy - 1, ix]:
|
||||
v[iz, iy, ix] = 0.0
|
||||
continue
|
||||
vc = v[iz, iy, ix]
|
||||
ve = v[iz, iy, ix + 1] if ix + 1 < nx else vc
|
||||
vw = v[iz, iy, ix - 1] if ix - 1 >= 0 else vc
|
||||
vn = v[iz, iy + 1, ix] if iy + 1 <= ny else vc
|
||||
vs = v[iz, iy - 1, ix]
|
||||
vt = v[iz + 1, iy, ix] if iz + 1 < nz else vc
|
||||
vb = v[iz - 1, iy, ix] if iz - 1 >= 0 else vc
|
||||
ue = 0.5 * (u[iz, iy, ix + 1] + u[iz, iy - 1, ix + 1]) if ix + 1 <= nx else 0.0
|
||||
uw = 0.5 * (u[iz, iy, ix] + u[iz, iy - 1, ix])
|
||||
wt = 0.5 * (w[iz + 1, iy, ix] + w[iz + 1, iy - 1, ix]) if iz + 1 < nz else 0.0
|
||||
wb = 0.5 * (w[iz, iy, ix] + w[iz, iy - 1, ix])
|
||||
dvy = (vc - vs) / dy if vc > 0 else (vn - vc) / dy
|
||||
dvx = (ve - vc) / dx if ue > 0 else (vc - vw) / dx
|
||||
dvz = (vt - vc) / dz if wt > 0 else (vc - vb) / dz
|
||||
conv = 0.5 * (ue + uw) * dvx + vc * dvy + 0.5 * (wt + wb) * dvz
|
||||
lap = (ve - 2.0 * vc + vw) / (dx * dx) + (vn - 2.0 * vc + vs) / (dy * dy) \
|
||||
+ (vt - 2.0 * vc + vb) / (dz * dz)
|
||||
vnew = vc + dt * (-conv + nu * lap)
|
||||
if plate[iz, iy, ix] or plate[iz, iy - 1, ix]:
|
||||
vnew /= (1.0 + dt * rporous)
|
||||
v[iz, iy, ix] = vnew
|
||||
# w faces: (nz+1, ny, nx)
|
||||
for iz in range(1, nz):
|
||||
for iy in range(ny):
|
||||
for ix in range(nx):
|
||||
if solid[iz, iy, ix] or solid[iz - 1, iy, ix]:
|
||||
w[iz, iy, ix] = 0.0
|
||||
continue
|
||||
wc = w[iz, iy, ix]
|
||||
we = w[iz, iy, ix + 1] if ix + 1 < nx else wc
|
||||
ww = w[iz, iy, ix - 1] if ix - 1 >= 0 else wc
|
||||
wn = w[iz, iy + 1, ix] if iy + 1 < ny else wc
|
||||
ws = w[iz, iy - 1, ix] if iy - 1 >= 0 else wc
|
||||
wt = w[iz + 1, iy, ix] if iz + 1 <= nz else wc
|
||||
wb = w[iz - 1, iy, ix]
|
||||
ue = 0.5 * (u[iz, iy, ix + 1] + u[iz - 1, iy, ix + 1]) if ix + 1 <= nx else 0.0
|
||||
uw = 0.5 * (u[iz, iy, ix] + u[iz - 1, iy, ix])
|
||||
vn = 0.5 * (v[iz, iy + 1, ix] + v[iz - 1, iy + 1, ix]) if iy + 1 <= ny else 0.0
|
||||
vs = 0.5 * (v[iz, iy, ix] + v[iz - 1, iy, ix])
|
||||
dwz = (wc - wb) / dz if wc > 0 else (wt - wc) / dz
|
||||
dwx = (we - wc) / dx if ue > 0 else (wc - ww) / dx
|
||||
dwy = (wn - wc) / dy if vn > 0 else (wc - ws) / dy
|
||||
conv = 0.5 * (ue + uw) * dwx + 0.5 * (vn + vs) * dwy + wc * dwz
|
||||
lap = (we - 2.0 * wc + ww) / (dx * dx) + (wn - 2.0 * wc + ws) / (dy * dy) \
|
||||
+ (wt - 2.0 * wc + wb) / (dz * dz)
|
||||
wnew = wc + dt * (-conv + nu * lap)
|
||||
if plate[iz, iy, ix] or plate[iz - 1, iy, ix]:
|
||||
wnew /= (1.0 + dt * rporous)
|
||||
w[iz, iy, ix] = wnew
|
||||
return
|
||||
|
||||
|
||||
@numba.njit(cache=True, fastmath=True)
|
||||
def apply_bc(u, v, w, solid, inlet, outlet, uin):
|
||||
nz, ny, nx = solid.shape
|
||||
# 压力进口: 边界面速度零梯度 (让气流自由进入)
|
||||
for iz in range(nz):
|
||||
for ix in range(nx):
|
||||
if inlet[iz, 0, ix]:
|
||||
v[iz, 0, ix] = v[iz, 1, ix]
|
||||
return
|
||||
|
||||
|
||||
@numba.njit(cache=True, fastmath=True)
|
||||
def divergence(u, v, w, solid, dx, dy, dz):
|
||||
nz, ny, nx = solid.shape
|
||||
d = np.zeros((nz, ny, nx))
|
||||
for iz in range(nz):
|
||||
for iy in range(ny):
|
||||
for ix in range(nx):
|
||||
if solid[iz, iy, ix]:
|
||||
continue
|
||||
d[iz, iy, ix] = (u[iz, iy, ix + 1] - u[iz, iy, ix]) / dx + \
|
||||
(v[iz, iy + 1, ix] - v[iz, iy, ix]) / dy + \
|
||||
(w[iz + 1, iy, ix] - w[iz, iy, ix]) / dz
|
||||
return d
|
||||
|
||||
|
||||
@numba.njit(cache=True, fastmath=True)
|
||||
def correct(u, v, w, p, solid, dt, rho, dx, dy, dz):
|
||||
nz, ny, nx = solid.shape
|
||||
c = dt / rho
|
||||
for iz in range(nz):
|
||||
for iy in range(ny):
|
||||
for ix in range(1, nx):
|
||||
if solid[iz, iy, ix] or solid[iz, iy, ix - 1]:
|
||||
continue
|
||||
u[iz, iy, ix] -= c * (p[iz, iy, ix] - p[iz, iy, ix - 1]) / dx
|
||||
for iz in range(nz):
|
||||
for iy in range(1, ny):
|
||||
for ix in range(nx):
|
||||
if solid[iz, iy, ix] or solid[iz, iy - 1, ix]:
|
||||
continue
|
||||
v[iz, iy, ix] -= c * (p[iz, iy, ix] - p[iz, iy - 1, ix]) / dy
|
||||
for iz in range(1, nz):
|
||||
for iy in range(ny):
|
||||
for ix in range(nx):
|
||||
if solid[iz, iy, ix] or solid[iz - 1, iy, ix]:
|
||||
continue
|
||||
w[iz, iy, ix] -= c * (p[iz, iy, ix] - p[iz - 1, iy, ix]) / dz
|
||||
return
|
||||
|
||||
|
||||
def build_poisson(solid, outlet, inlet):
|
||||
n = NX * NY * NZ
|
||||
idx = np.arange(n).reshape(NZ, NY, NX)
|
||||
fluid = ~solid
|
||||
rows, cols, vals = [], [], []
|
||||
|
||||
def add(r, c, v):
|
||||
rows.append(np.asarray(r).ravel()); cols.append(np.asarray(c).ravel())
|
||||
vals.append(np.asarray(v).ravel())
|
||||
|
||||
diag = np.zeros((NZ, NY, NX))
|
||||
m = fluid[:, :, :-1] & fluid[:, :, 1:]
|
||||
c = 1.0 / DX ** 2
|
||||
add(idx[:, :, :-1][m], idx[:, :, 1:][m], np.full(m.sum(), -c))
|
||||
add(idx[:, :, 1:][m], idx[:, :, :-1][m], np.full(m.sum(), -c))
|
||||
diag[:, :, :-1] += m * c; diag[:, :, 1:] += m * c
|
||||
m = fluid[:, :-1, :] & fluid[:, 1:, :]
|
||||
c = 1.0 / DY ** 2
|
||||
add(idx[:, :-1, :][m], idx[:, 1:, :][m], np.full(m.sum(), -c))
|
||||
add(idx[:, 1:, :][m], idx[:, :-1, :][m], np.full(m.sum(), -c))
|
||||
diag[:, :-1, :] += m * c; diag[:, 1:, :] += m * c
|
||||
m = fluid[:-1, :, :] & fluid[1:, :, :]
|
||||
c = 1.0 / DZ ** 2
|
||||
add(idx[:-1, :, :][m], idx[1:, :, :][m], np.full(m.sum(), -c))
|
||||
add(idx[1:, :, :][m], idx[:-1, :, :][m], np.full(m.sum(), -c))
|
||||
diag[:-1, :, :] += m * c; diag[1:, :, :] += m * c
|
||||
|
||||
fl = fluid.ravel()
|
||||
df = diag.ravel().copy()
|
||||
df = np.where(fl & (df == 0.0), 1.0, df)
|
||||
df[~fl] = 1.0
|
||||
add(np.arange(n), np.arange(n), df)
|
||||
A = csr_matrix((np.concatenate(vals), (np.concatenate(rows), np.concatenate(cols))),
|
||||
shape=(n, n))
|
||||
of = outlet.ravel()
|
||||
inf = inlet.ravel()
|
||||
fix = of | inf
|
||||
if fix.any():
|
||||
oi = np.where(fix)[0]
|
||||
A = A.tolil(); A[oi, :] = 0.0; A[oi, oi] = 1.0; A = A.tocsr()
|
||||
return A, idx
|
||||
|
||||
|
||||
def main():
|
||||
t_all = time.time()
|
||||
print(f"[rans] staggered grid {NX}x{NY}x{NZ}={NX*NY*NZ} PPM={PPM}")
|
||||
solid, plate, inlet, outlet = build_masks()
|
||||
print(f"[rans] fluid={int((~solid).sum())} plate={int(plate.sum())} R_porous={R_POROUS:.0f}")
|
||||
|
||||
A, idx = build_poisson(solid, outlet, inlet)
|
||||
print(f"[rans] poisson nnz={A.nnz}, AMG ...")
|
||||
t0 = time.time()
|
||||
ml = pyamg.smoothed_aggregation_solver(A, max_coarse=300)
|
||||
print(f"[rans] AMG {time.time()-t0:.2f}s levels={len(ml.levels)}")
|
||||
|
||||
u = np.zeros((NZ, NY, NX + 1))
|
||||
v = np.zeros((NZ, NY + 1, NX))
|
||||
w = np.zeros((NZ + 1, NY, NX))
|
||||
p = np.zeros((NZ, NY, NX))
|
||||
dt = CFL * min(DX, DY, DZ) / (U_IN + 1e-9)
|
||||
print(f"[rans] 恒功率风机: P={P_FAN:.0f} W 额定({Q_RATED*3600:.0f} m3/h, {DP_RATED:.0f} Pa) "
|
||||
f"关死={DP_SHUT:.0f} Pa nu_eff={NU_EFF:.4f} dt={dt:.3e}")
|
||||
|
||||
p_fan = DP_RATED
|
||||
t0 = time.time()
|
||||
for it in range(MAX_STEPS):
|
||||
predict(u, v, w, solid, plate, inlet, outlet, NU_EFF, dt, DX, DY, DZ, R_POROUS, U_IN)
|
||||
d = divergence(u, v, w, solid, DX, DY, DZ)
|
||||
rhs = (RHO / dt) * d
|
||||
rhs[outlet] = 0.0
|
||||
rhs[inlet] = p_fan
|
||||
sol = ml.solve(rhs.ravel(), x0=p.ravel(), tol=1e-8, maxiter=200)
|
||||
p = sol.reshape(NZ, NY, NX)
|
||||
correct(u, v, w, p, solid, dt, RHO, DX, DY, DZ)
|
||||
apply_bc(u, v, w, solid, inlet, outlet, U_IN)
|
||||
# 测进口流量 -> 更新风机压力
|
||||
iz_in = int(round((-0.45 - Y0) / DY))
|
||||
m_in = ~solid[:, iz_in, :]
|
||||
q_in = abs(v[:, iz_in, :][m_in].sum()) * DX * DZ
|
||||
p_fan = fan_dp(q_in)
|
||||
if it % 100 == 0:
|
||||
uc = 0.5 * (u[:, :, :-1] + u[:, :, 1:])
|
||||
vc = 0.5 * (v[:, :-1, :] + v[:, 1:, :])
|
||||
wc = 0.5 * (w[:-1, :, :] + w[1:, :, :])
|
||||
sp = np.sqrt(uc ** 2 + vc ** 2 + wc ** 2)
|
||||
print(f" it={it:5d} t={time.time()-t0:.1f}s |U|max={sp.max():.2f} "
|
||||
f"Q_in={q_in*3600:.0f} p_fan={p_fan:.0f} p=[{p.min():.0f},{p.max():.0f}]")
|
||||
if time.time() - t0 > TIME_BUDGET:
|
||||
break
|
||||
|
||||
izp = int(round((PLATE_T - Z0) / DZ))
|
||||
q_plate = w[izp][~solid[izp]].sum() * DX * DY * 3600.0
|
||||
iz_in = int(round((-0.45 - Y0) / DY))
|
||||
m_in = ~solid[:, iz_in, :]
|
||||
q_in = v[:, iz_in, :][m_in].sum() * DX * DZ * 3600.0
|
||||
iz_out = NZ - 2
|
||||
q_out = w[iz_out][~solid[iz_out]].sum() * DX * DY * 3600.0
|
||||
a0, a1 = int(round((-0.30 - Z0) / DZ)), int(round((-0.05 - Z0) / DZ))
|
||||
p_pl = p[a0:a1][~solid[a0:a1]].mean()
|
||||
b0, b1 = int(round((0.15 - Z0) / DZ)), int(round((0.35 - Z0) / DZ))
|
||||
p_fb = p[b0:b1][~solid[b0:b1]].mean()
|
||||
print(f"[rans] steps={it+1} solve={time.time()-t0:.1f}s")
|
||||
print(f"[rans] plate Q={q_plate:.0f} inlet Q={q_in:.0f} outlet Q={q_out:.0f} m3/h")
|
||||
print(f"[rans] p_风箱={p_pl:.0f} p_沸腾段={p_fb:.0f} Δp={p_pl-p_fb:+.0f} Pa")
|
||||
out = os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))), "docs", "assets")
|
||||
np.savez_compressed(os.path.join(out, "fbed_rans.npz"), u=u, v=v, w=w, p=p,
|
||||
solid=solid, plate=plate)
|
||||
print(f"[rans] total {time.time()-t_all:.1f}s")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -0,0 +1,507 @@
|
||||
#!/usr/bin/env python
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
水平卧式流化床 —— 空床气相流场快速仿真 (3D D3Q19 LBM, GPU)
|
||||
============================================================
|
||||
几何 (m)
|
||||
--------
|
||||
布风板 : 2400 (x, 长度) x 567 (y, 宽度), 水平
|
||||
孔板 : 孔径 2 mm, 开孔率 5.6%, 孔轴向 +x 倾斜 5 度 (朝出料端)
|
||||
进风分区 : 沿长度均分 4 段 (4 x 600 mm), 每段下方独立风箱
|
||||
风箱 : 每段 600 x 567 x 450 mm, DN300 进风管从 567 侧面水平居中进入
|
||||
沸腾/扩大段 : 由 2400 x 567 在 1000 mm 高度内矩形渐扩至 2400 x 1200
|
||||
上部稳压段 : 2400 x 1200 x 400 mm
|
||||
排风 : 顶部居中 DN800 圆管
|
||||
|
||||
工况
|
||||
----
|
||||
单段风量 : 2500 m^3/h (4 段共 10000 m^3/h)
|
||||
床板压差 : 2900 Pa (布风板上下)
|
||||
表观气速 : U_sup = 10000/3600 / (2.4*0.567) = 2.04 m/s
|
||||
管口气速 : DN300 单管 9.82 m/s ; 排风 DN800 5.53 m/s
|
||||
|
||||
方法
|
||||
----
|
||||
D3Q19 LBM (局部显式 stream+collide, GPU/numba.cuda), 目标 < 20 s。
|
||||
布风板用"多孔执行面"模型: 以表观气速 (含 5 度倾斜分量) 均匀喷射, 质量守恒。
|
||||
进风管 / 排风管用速度边界 (velocity inlet / outlet)。
|
||||
2900 Pa 远超声速不可压 LBM 可表达范围, 故用"给定风量"而非"给定压差"驱动。
|
||||
|
||||
输出
|
||||
----
|
||||
docs/assets/fbed-overview.png 四联总览
|
||||
docs/assets/fbed-slice-xz.png 垂直切面 y=0 (风箱-布风板-沸腾段-排风)
|
||||
docs/assets/fbed-slice-yz.png 垂直切面 x=1.2 (侧向渐扩)
|
||||
docs/assets/fbed-slice-xy-plate.png 水平切面 z=0.15 (布风板近区射流)
|
||||
docs/assets/fbed-slice-xy-free.png 水平切面 z=0.9 (沸腾段)
|
||||
|
||||
运行
|
||||
----
|
||||
python tools/fbed_simulation.py
|
||||
"""
|
||||
import os
|
||||
import time
|
||||
|
||||
import numpy as np
|
||||
from numba import cuda
|
||||
|
||||
import matplotlib
|
||||
matplotlib.use("Agg")
|
||||
import matplotlib.pyplot as plt
|
||||
from matplotlib.colors import LinearSegmentedColormap
|
||||
|
||||
# --------------------------------------------------------------------------
|
||||
# 几何参数 (m)
|
||||
# --------------------------------------------------------------------------
|
||||
BED_L, BED_W = 2.4, 0.567
|
||||
N_SEC = 4
|
||||
SEC_L = BED_L / N_SEC
|
||||
PLENUM_H = 0.45
|
||||
PLATE_T = 0.002
|
||||
|
||||
FB_H, FB_W = 1.0, 1.2
|
||||
UPPER_H = 0.4
|
||||
REDUCER_H = 0.5
|
||||
DUCT_H = 0.3
|
||||
EXH_R = 0.4
|
||||
|
||||
PIPE_R = 0.15
|
||||
PIPE_Y0 = -0.65
|
||||
|
||||
X0, X1 = 0.0, BED_L
|
||||
Y0, Y1 = PIPE_Y0, 0.65
|
||||
Z1 = FB_H + UPPER_H + REDUCER_H + DUCT_H
|
||||
Z0 = -PLENUM_H
|
||||
PPM = int(os.environ.get("FBED_PPM", "100"))
|
||||
NX = int(round((X1 - X0) * PPM))
|
||||
NY = int(round((Y1 - Y0) * PPM))
|
||||
NZ = int(round((Z1 - Z0) * PPM))
|
||||
DX = 1.0 / PPM
|
||||
|
||||
# --------------------------------------------------------------------------
|
||||
# 工况
|
||||
# --------------------------------------------------------------------------
|
||||
Q_SEC = 2500.0 / 3600.0
|
||||
A_BED = BED_L * BED_W
|
||||
U_SUP = Q_SEC * N_SEC / A_BED
|
||||
TILT = np.deg2rad(5.0)
|
||||
TILT_DIR = -1.0 # 绿色箭头方向: 孔轴向上-左 (-x)
|
||||
UP_X = TILT_DIR * U_SUP * np.sin(TILT)
|
||||
UP_Z = U_SUP * np.cos(TILT)
|
||||
TILT_TAN = np.tan(TILT)
|
||||
A_PIPE = np.pi * PIPE_R ** 2
|
||||
U_IN = Q_SEC / A_PIPE
|
||||
A_EXH = np.pi * EXH_R ** 2
|
||||
U_OUT = Q_SEC * N_SEC / A_EXH
|
||||
|
||||
U_LB = float(os.environ.get("FBED_ULB", "0.03"))
|
||||
SCALE_U = U_IN / U_LB
|
||||
U_IN_LB = U_LB
|
||||
U_OUT_LB = U_OUT / SCALE_U
|
||||
UPX_LB = UP_X / SCALE_U
|
||||
UPZ_LB = UP_Z / SCALE_U
|
||||
NU_LB = 0.04
|
||||
TAU = 0.5 + 3.0 * NU_LB
|
||||
|
||||
DP_RATED = float(os.environ.get("FBED_DP", "2900.0"))
|
||||
Q_RATED = (2500.0 * 4) / 3600.0
|
||||
P_FAN = DP_RATED * Q_RATED
|
||||
DP_SHUT = DP_RATED * 1.5
|
||||
|
||||
|
||||
def fan_dp(q_m3s):
|
||||
"""恒功率风机曲线 dp*Q=P -> dp=P/Q (关死压力截断)"""
|
||||
if q_m3s < 1e-6:
|
||||
return DP_SHUT
|
||||
return min(P_FAN / q_m3s, DP_SHUT)
|
||||
|
||||
TIME_BUDGET = float(os.environ.get("FBED_TB", "120.0"))
|
||||
MAX_STEPS = 40000
|
||||
|
||||
EX = np.array([0, 1, -1, 0, 0, 0, 0, 1, -1, 1, -1, 1, -1, 1, -1, 0, 0, 0, 0], dtype=np.int32)
|
||||
EY = np.array([0, 0, 0, 1, -1, 0, 0, 1, -1, -1, 1, 0, 0, 0, 0, 1, -1, 1, -1], dtype=np.int32)
|
||||
EZ = np.array([0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 1, -1, -1, 1, 1, -1, -1, 1], dtype=np.int32)
|
||||
WT = np.array([1/3] + [1/18]*6 + [1/36]*12, dtype=np.float64)
|
||||
OPP = np.array([0, 2, 1, 4, 3, 6, 5, 8, 7, 10, 9, 12, 11, 14, 13, 16, 15, 18, 17], dtype=np.int32)
|
||||
|
||||
XC = X0 + (np.arange(NX) + 0.5) / PPM
|
||||
YC = Y0 + (np.arange(NY) + 0.5) / PPM
|
||||
ZC = Z0 + (np.arange(NZ) + 0.5) / PPM
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------
|
||||
# GPU 核函数
|
||||
# --------------------------------------------------------------------------
|
||||
@cuda.jit
|
||||
def collide_kernel(f, rho, ux, uy, uz, solid, plate, inlet, outlet,
|
||||
ex, ey, ez, w, tau, drho, ttilt, nz, ny, nx):
|
||||
idx = cuda.grid(1)
|
||||
if idx >= nz * ny * nx:
|
||||
return
|
||||
ix = idx % nx
|
||||
iy = (idx // nx) % ny
|
||||
iz = idx // (nx * ny)
|
||||
if solid[iz, iy, ix]:
|
||||
return
|
||||
r = 0.0
|
||||
mx = 0.0
|
||||
my = 0.0
|
||||
mz = 0.0
|
||||
for k in range(19):
|
||||
fk = f[k, iz, iy, ix]
|
||||
r += fk
|
||||
mx += ex[k] * fk
|
||||
my += ey[k] * fk
|
||||
mz += ez[k] * fk
|
||||
inv = 1.0 / r
|
||||
u = mx * inv
|
||||
v = my * inv
|
||||
ww = mz * inv
|
||||
om = 1.0 / tau
|
||||
uu = u * u + v * v + ww * ww
|
||||
for k in range(19):
|
||||
eu = ex[k] * u + ey[k] * v + ez[k] * ww
|
||||
feq = w[k] * r * (1.0 + 3.0 * eu + 4.5 * eu * eu - 1.5 * uu)
|
||||
f[k, iz, iy, ix] -= (f[k, iz, iy, ix] - feq) * om
|
||||
if plate[iz, iy, ix]:
|
||||
# 布风板 5 度倾斜: 开孔处施加水平动量 (绿箭头方向 -x)
|
||||
fx = -ttilt * r * ww
|
||||
for k in range(19):
|
||||
f[k, iz, iy, ix] += 3.0 * w[k] * ex[k] * fx
|
||||
u += fx * inv
|
||||
if inlet[iz, iy, ix]:
|
||||
# 压力进口: 固定密度(= 风机压力), 速度零梯度自由
|
||||
r = 1.0 + drho
|
||||
uu = u * u + v * v + ww * ww
|
||||
for k in range(19):
|
||||
eu = ex[k] * u + ey[k] * v + ez[k] * ww
|
||||
f[k, iz, iy, ix] = w[k] * r * (1.0 + 3.0 * eu + 4.5 * eu * eu - 1.5 * uu)
|
||||
elif outlet[iz, iy, ix]:
|
||||
r = 1.0
|
||||
uu = u * u + v * v + ww * ww
|
||||
for k in range(19):
|
||||
eu = ex[k] * u + ey[k] * v + ez[k] * ww
|
||||
f[k, iz, iy, ix] = w[k] * r * (1.0 + 3.0 * eu + 4.5 * eu * eu - 1.5 * uu)
|
||||
rho[iz, iy, ix] = r
|
||||
ux[iz, iy, ix] = u
|
||||
uy[iz, iy, ix] = v
|
||||
uz[iz, iy, ix] = ww
|
||||
|
||||
|
||||
@cuda.jit
|
||||
def stream_kernel(f, ftmp, solid, ex, ey, ez, opp, nz, ny, nx):
|
||||
idx = cuda.grid(1)
|
||||
if idx >= nz * ny * nx:
|
||||
return
|
||||
ix = idx % nx
|
||||
iy = (idx // nx) % ny
|
||||
iz = idx // (nx * ny)
|
||||
if solid[iz, iy, ix]:
|
||||
return
|
||||
for k in range(19):
|
||||
sz = iz - ez[k]
|
||||
sy = iy - ey[k]
|
||||
sx = ix - ex[k]
|
||||
if 0 <= sz < nz and 0 <= sy < ny and 0 <= sx < nx and not solid[sz, sy, sx]:
|
||||
ftmp[k, iz, iy, ix] = f[k, sz, sy, sx]
|
||||
else:
|
||||
ftmp[k, iz, iy, ix] = f[opp[k], iz, iy, ix]
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------
|
||||
# 几何
|
||||
# --------------------------------------------------------------------------
|
||||
def _xi(v):
|
||||
return int(round((v - X0) * PPM))
|
||||
|
||||
|
||||
def _yi(v):
|
||||
return int(round((v - Y0) * PPM))
|
||||
|
||||
|
||||
def _zi(v):
|
||||
return int(round((v - Z0) * PPM))
|
||||
|
||||
|
||||
def _clip(i, n):
|
||||
return max(0, min(i, n))
|
||||
|
||||
|
||||
def carve_box(solid, xa, xb, ya, yb, za, zb):
|
||||
ix0, ix1 = _clip(_xi(xa), NX), _clip(_xi(xb), NX)
|
||||
iy0, iy1 = _clip(_yi(ya), NY), _clip(_yi(yb), NY)
|
||||
iz0, iz1 = _clip(_zi(za), NZ), _clip(_zi(zb), NZ)
|
||||
solid[iz0:iz1, iy0:iy1, ix0:ix1] = False
|
||||
|
||||
|
||||
def build_geometry():
|
||||
solid = np.ones((NZ, NY, NX), dtype=np.bool_)
|
||||
yw = BED_W / 2.0
|
||||
|
||||
carve_box(solid, 0.0, BED_L, -yw, yw, -PLENUM_H, 0.0)
|
||||
for k in range(1, N_SEC):
|
||||
xk = k * SEC_L
|
||||
carve_box(solid, xk - 0.0125, xk + 0.0125, -yw, yw, -PLENUM_H, 0.0)
|
||||
for k in range(1, N_SEC):
|
||||
xk = k * SEC_L
|
||||
solid[_zi(-PLENUM_H):_zi(0.0), _yi(-yw):_yi(yw), _xi(xk - 0.0125):_xi(xk + 0.0125)] = True
|
||||
|
||||
carve_box(solid, 0.0, BED_L, -yw, yw, 0.0, PLATE_T)
|
||||
|
||||
for iz in range(_zi(PLATE_T), _zi(FB_H)):
|
||||
z = ZC[iz]
|
||||
w = yw + (FB_W / 2.0 - yw) * (z / FB_H)
|
||||
solid[iz, :, :] = np.where(
|
||||
((YC[:, None] >= -w) & (YC[:, None] <= w) &
|
||||
(XC[None, :] >= 0.0) & (XC[None, :] <= BED_L)),
|
||||
False, solid[iz, :, :])
|
||||
carve_box(solid, 0.0, BED_L, -FB_W / 2.0, FB_W / 2.0, FB_H, FB_H + UPPER_H)
|
||||
|
||||
z_r0 = FB_H + UPPER_H
|
||||
z_r1 = z_r0 + REDUCER_H
|
||||
for iz in range(_zi(z_r0), _zi(z_r1)):
|
||||
t = (ZC[iz] - z_r0) / REDUCER_H
|
||||
xh = BED_L / 2.0 - (BED_L / 2.0 - 0.4) * t
|
||||
yh = FB_W / 2.0 - (FB_W / 2.0 - 0.4) * t
|
||||
solid[iz] = np.where(
|
||||
(np.abs(XC[None, :] - BED_L / 2.0) <= xh) & (np.abs(YC[:, None]) <= yh),
|
||||
False, solid[iz])
|
||||
|
||||
cyl = (XC[None, :] - BED_L / 2.0) ** 2 + (YC[:, None] - 0.0) ** 2 <= EXH_R ** 2
|
||||
iz_e0, iz_e1 = _zi(z_r1), NZ
|
||||
solid[iz_e0:iz_e1] = np.where(cyl[None, :, :], False, solid[iz_e0:iz_e1])
|
||||
|
||||
iy_p0, iy_p1 = _yi(PIPE_Y0), _yi(-yw)
|
||||
for k in range(N_SEC):
|
||||
xk = SEC_L * (k + 0.5)
|
||||
zk = -PLENUM_H / 2.0
|
||||
m = ((XC[None, :] - xk) ** 2 + (ZC[:, None] - zk) ** 2) <= PIPE_R ** 2
|
||||
m3 = np.repeat(m[:, None, :], iy_p1 - iy_p0, axis=1)
|
||||
region = solid[:, iy_p0:iy_p1, :]
|
||||
region[m3] = False
|
||||
solid[:, iy_p0:iy_p1, :] = region
|
||||
|
||||
iz_plate = _zi(0.0)
|
||||
# 多孔布风板: 固体 + 规则开孔 (孔 hs x hs 格, 间距 ss 格)
|
||||
izp0, izp1 = _zi(0.0), max(_zi(PLATE_T), _zi(0.0) + 1)
|
||||
iy0p, iy1p = _yi(-yw), _yi(yw)
|
||||
solid[izp0:izp1, iy0p:iy1p, :] = True
|
||||
hs = int(os.environ.get("FBED_HOLE", "1"))
|
||||
ss = int(os.environ.get("FBED_SPACE", "4"))
|
||||
for iz in range(izp0, izp1):
|
||||
for iy in range(iy0p, iy1p):
|
||||
for ix in range(NX):
|
||||
if ((iy - iy0p) % ss < hs) and (ix % ss < hs):
|
||||
solid[iz, iy, ix] = False
|
||||
plate = np.zeros_like(solid)
|
||||
plate[izp0:izp1, iy0p:iy1p, :] = ~solid[izp0:izp1, iy0p:iy1p, :]
|
||||
|
||||
inlet = np.zeros_like(solid)
|
||||
inlet[:, iy_p0:iy_p0 + 3, :] = ~solid[:, iy_p0:iy_p0 + 3, :]
|
||||
|
||||
outlet = np.zeros_like(solid)
|
||||
outlet[NZ - 2:, :, :] = ~solid[NZ - 2:, :, :]
|
||||
return solid, plate, inlet, outlet
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------
|
||||
# 求解
|
||||
# --------------------------------------------------------------------------
|
||||
def run_solver(solid, plate, inlet, outlet):
|
||||
f = np.zeros((19, NZ, NY, NX), dtype=np.float32)
|
||||
for k in range(19):
|
||||
f[k] = WT[k]
|
||||
ftmp = np.zeros_like(f)
|
||||
rho = np.ones((NZ, NY, NX), dtype=np.float32)
|
||||
ux = np.zeros_like(rho)
|
||||
uy = np.zeros_like(rho)
|
||||
uz = np.zeros_like(rho)
|
||||
|
||||
d = {name: cuda.to_device(a) for name, a in [
|
||||
("f", f), ("ftmp", ftmp), ("rho", rho), ("ux", ux), ("uy", uy), ("uz", uz),
|
||||
("solid", solid), ("plate", plate), ("inlet", inlet), ("outlet", outlet),
|
||||
("ex", EX), ("ey", EY), ("ez", EZ), ("w", WT), ("opp", OPP)]}
|
||||
|
||||
n = NX * NY * NZ
|
||||
tpb = 128
|
||||
bpg = (n + tpb - 1) // tpb
|
||||
|
||||
drho_fan = 3.0 * DP_RATED / (1.2 * SCALE_U ** 2)
|
||||
cuda.synchronize()
|
||||
t0 = time.time()
|
||||
step = 0
|
||||
for step in range(MAX_STEPS):
|
||||
collide_kernel[bpg, tpb](d["f"], d["rho"], d["ux"], d["uy"], d["uz"],
|
||||
d["solid"], d["plate"], d["inlet"], d["outlet"],
|
||||
d["ex"], d["ey"], d["ez"], d["w"], TAU,
|
||||
drho_fan, 0.0, NZ, NY, NX)
|
||||
stream_kernel[bpg, tpb](d["f"], d["ftmp"], d["solid"],
|
||||
d["ex"], d["ey"], d["ez"], d["opp"], NZ, NY, NX)
|
||||
d["f"], d["ftmp"] = d["ftmp"], d["f"]
|
||||
if step % 100 == 99:
|
||||
cuda.synchronize()
|
||||
vy = d["uy"].copy_to_host()
|
||||
iz_in = _yi(-0.45)
|
||||
q_in = float(vy[:, iz_in, :][~solid[:, iz_in, :]].sum()) * SCALE_U * DX * DX
|
||||
dp = fan_dp(abs(q_in))
|
||||
drho_fan = 3.0 * dp / (1.2 * SCALE_U ** 2)
|
||||
if step % 2000 == 1999:
|
||||
print(f" step {step+1:6d} Q_in={abs(q_in)*3600:.0f} m3/h "
|
||||
f"dp_fan={dp:.0f} Pa t={time.time()-t0:.1f}s")
|
||||
if time.time() - t0 > TIME_BUDGET:
|
||||
break
|
||||
cuda.synchronize()
|
||||
solve_dt = time.time() - t0
|
||||
steps = step + 1
|
||||
return (d["ux"].copy_to_host(), d["uy"].copy_to_host(), d["uz"].copy_to_host(),
|
||||
d["rho"].copy_to_host(), steps, solve_dt)
|
||||
|
||||
|
||||
def metrics(ux, uy, uz, rho, solid, plate, inlet, outlet, steps, solve_dt):
|
||||
fluid = ~solid
|
||||
speed = np.sqrt(ux ** 2 + uy ** 2 + uz ** 2) * SCALE_U
|
||||
umax = float(speed[fluid].max())
|
||||
iz_max, iy_max, ix_max = np.unravel_index(np.argmax(np.where(fluid, speed, 0)), speed.shape)
|
||||
|
||||
cell_a = DX * DX
|
||||
iz_plate = _zi(PLATE_T)
|
||||
fl_p = ~solid[iz_plate]
|
||||
q_plate = float(uz[iz_plate][fl_p].sum()) * SCALE_U * cell_a * 3600.0
|
||||
iy_in = _yi(-0.45)
|
||||
q_in = float(uy[:, iy_in, :][~solid[:, iy_in, :]].sum()) * SCALE_U * cell_a * 3600.0
|
||||
q_out = float(uz[NZ - 2][outlet[NZ - 2]].sum()) * SCALE_U * cell_a * 3600.0
|
||||
|
||||
p = (rho - 1.0) / 3.0 * (1.2 * SCALE_U ** 2)
|
||||
m_pl = ~solid[_zi(-0.30):_zi(-0.05)]
|
||||
p_pl = p[_zi(-0.30):_zi(-0.05)][m_pl].mean()
|
||||
m_fb = ~solid[_zi(PLATE_T + 0.05):_zi(PLATE_T + 0.25)]
|
||||
p_fb = p[_zi(PLATE_T + 0.05):_zi(PLATE_T + 0.25)][m_fb].mean()
|
||||
r_fluid = rho[~solid]
|
||||
print(f"[fbed] grid {NX}x{NY}x{NZ}={NX*NY*NZ}, steps={steps}, solve={solve_dt:.2f} s")
|
||||
print(f"[fbed] |U|max={umax:.2f} m/s @ (x={XC[ix_max]:.2f}, y={YC[iy_max]:.2f}, z={ZC[iz_max]:.2f})")
|
||||
print(f"[fbed] plate Q={q_plate:.0f} inlet Q={q_in:.0f} outlet Q={q_out:.0f} m3/h (target 10000)")
|
||||
print(f"[fbed] p_风箱={p_pl:.0f} Pa p_沸腾段={p_fb:.0f} Pa Δp(风箱-沸腾段)={p_pl-p_fb:+.0f} Pa")
|
||||
print(f"[fbed] rho=[{r_fluid.min():.4f},{r_fluid.max():.4f}] p99={np.percentile(r_fluid,99):.4f}")
|
||||
return speed, p
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------
|
||||
# 绘图
|
||||
# --------------------------------------------------------------------------
|
||||
BG = "#0d1117"
|
||||
FG = "#e6edf3"
|
||||
MUTED = "#8b949e"
|
||||
ACCENT = "#58a6ff"
|
||||
CMAP_FLOW = LinearSegmentedColormap.from_list(
|
||||
"flow", ["#0a1626", "#123c63", "#1f6feb", "#39d0d8", "#7ee787", "#f0f6fc"])
|
||||
plt.rcParams.update({
|
||||
"font.sans-serif": ["Microsoft YaHei", "SimHei", "DejaVu Sans"],
|
||||
"axes.unicode_minus": False,
|
||||
"figure.facecolor": BG, "savefig.facecolor": BG, "axes.facecolor": BG,
|
||||
"text.color": FG, "axes.labelcolor": FG, "xtick.color": MUTED,
|
||||
"ytick.color": MUTED, "axes.edgecolor": "#30363d",
|
||||
})
|
||||
|
||||
OUT = os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))),
|
||||
"docs", "assets")
|
||||
|
||||
|
||||
def plane(speed, ux, uy, uz, solid, which, idx):
|
||||
if which == "xz":
|
||||
return (speed[:, idx, :], ux[:, idx, :], uz[:, idx, :], solid[:, idx, :],
|
||||
[X0, X1, Z0, Z1], "x (m)", "z (m)")
|
||||
if which == "yz":
|
||||
return (speed[:, :, idx], uy[:, :, idx], uz[:, :, idx], solid[:, :, idx],
|
||||
[Y0, Y1, Z0, Z1], "y (m)", "z (m)")
|
||||
return (speed[idx, :, :], ux[idx, :, :], uy[idx, :, :], solid[idx, :, :],
|
||||
[X0, X1, Y0, Y1], "x (m)", "y (m)")
|
||||
|
||||
|
||||
def draw_plane(ax, data, vmax, density=1.1):
|
||||
s2, u2, v2, sol2, ext, xl, yl = data
|
||||
masked = np.where(sol2, np.nan, s2)
|
||||
im = ax.imshow(masked, origin="lower", extent=ext, cmap=CMAP_FLOW,
|
||||
vmin=0, vmax=vmax, interpolation="bilinear", aspect="equal")
|
||||
ny, nx = sol2.shape
|
||||
xs = np.linspace(ext[0], ext[1], nx)
|
||||
ys = np.linspace(ext[2], ext[3], ny)
|
||||
U = np.where(sol2, np.nan, u2)
|
||||
V = np.where(sol2, np.nan, v2)
|
||||
with np.errstate(invalid="ignore"):
|
||||
ax.streamplot(xs, ys, U, V, density=density, color=(0.85, 0.88, 0.92, 0.4),
|
||||
linewidth=0.5, arrowsize=0.6, zorder=5)
|
||||
ax.contour(np.linspace(ext[0], ext[1], sol2.shape[1]),
|
||||
np.linspace(ext[2], ext[3], sol2.shape[0]),
|
||||
sol2.astype(float), levels=[0.5], colors="#adbac7", linewidths=0.8, zorder=6)
|
||||
ax.set_xlabel(xl, fontsize=10)
|
||||
ax.set_ylabel(yl, fontsize=10)
|
||||
ax.set_xlim(ext[0], ext[1])
|
||||
ax.set_ylim(ext[2], ext[3])
|
||||
return im
|
||||
|
||||
|
||||
def save_slice(speed, ux, uy, uz, solid, which, idx, title, name, vmax):
|
||||
data = plane(speed, ux, uy, uz, solid, which, idx)
|
||||
fig, ax = plt.subplots(figsize=(10.5, 7.6))
|
||||
im = draw_plane(ax, data, vmax)
|
||||
cb = fig.colorbar(im, ax=ax, fraction=0.036, pad=0.02)
|
||||
cb.set_label("速度大小 |U| (m/s)", fontsize=10)
|
||||
cb.outline.set_edgecolor("#30363d")
|
||||
ax.set_title(title, fontsize=13, pad=10)
|
||||
fig.tight_layout()
|
||||
fig.savefig(os.path.join(OUT, name), dpi=540)
|
||||
plt.close(fig)
|
||||
print(" ->", name)
|
||||
|
||||
|
||||
def save_overview(speed, ux, uy, uz, solid, steps, solve_dt, umax):
|
||||
fig, axes = plt.subplots(2, 2, figsize=(16.0, 12.0))
|
||||
specs = [
|
||||
("xz", _yi(0.0), "(a) 垂直切面 y=0 —— 风箱 / 布风板 / 沸腾段 / 排风", 4.0),
|
||||
("yz", _xi(BED_L / 2.0), "(b) 垂直切面 x=1.2 —— 侧向渐扩（567→1200 mm）", 4.0),
|
||||
("xy", _zi(0.15), "(c) 水平切面 z=0.15 —— 布风板近区射流", 3.0),
|
||||
("xy", _zi(0.90), "(d) 水平切面 z=0.90 —— 沸腾段", 3.0),
|
||||
]
|
||||
for ax, (which, idx, title, vmax) in zip(axes.ravel(), specs):
|
||||
data = plane(speed, ux, uy, uz, solid, which, idx)
|
||||
im = draw_plane(ax, data, vmax)
|
||||
cb = fig.colorbar(im, ax=ax, fraction=0.036, pad=0.02)
|
||||
cb.set_label("|U| (m/s)", fontsize=9.5)
|
||||
cb.outline.set_edgecolor("#30363d")
|
||||
cb.ax.tick_params(labelsize=8.5)
|
||||
ax.set_title(title, fontsize=12, pad=8)
|
||||
fig.suptitle("水平卧式流化床 空床流场 —— 4 段侧进风 / 5.6% 开孔倾斜布风板 / 顶部 DN800 排风",
|
||||
fontsize=16, y=0.985)
|
||||
fig.text(0.5, 0.006,
|
||||
f"求解 {solve_dt:.1f} s · {steps} 步 · |U|max={umax:.1f} m/s · "
|
||||
f"总风量 10000 m³/h · D3Q19 LBM (GPU)",
|
||||
ha="center", fontsize=10.5, color=MUTED)
|
||||
fig.tight_layout(rect=[0, 0.02, 1, 0.97])
|
||||
fig.savefig(os.path.join(OUT, "fbed-overview.png"), dpi=520)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-overview.png")
|
||||
|
||||
|
||||
def main():
|
||||
os.makedirs(OUT, exist_ok=True)
|
||||
print(f"[fbed] U_sup={U_SUP:.3f} m/s U_in(DN300)={U_IN:.2f} m/s "
|
||||
f"U_out(DN800)={U_OUT:.2f} m/s tilt=5 deg")
|
||||
print(f"[fbed] 恒功率风机: P={P_FAN:.0f} W 额定({Q_RATED*3600:.0f} m3/h, {DP_RATED:.0f} Pa) "
|
||||
f"关死={DP_SHUT:.0f} Pa")
|
||||
solid, plate, inlet, outlet = build_geometry()
|
||||
ux, uy, uz, rho, steps, solve_dt = run_solver(solid, plate, inlet, outlet)
|
||||
speed, p = metrics(ux, uy, uz, rho, solid, plate, inlet, outlet, steps, solve_dt)
|
||||
|
||||
np.savez_compressed(os.path.join(OUT, "fbed_fields.npz"),
|
||||
ux=ux, uy=uy, uz=uz, rho=rho, solid=solid,
|
||||
plate=plate, inlet=inlet, outlet=outlet,
|
||||
steps=steps, solve_dt=solve_dt)
|
||||
print(" -> fbed_fields.npz")
|
||||
|
||||
t0 = time.time()
|
||||
save_overview(speed, ux, uy, uz, solid, steps, solve_dt, float(speed[~solid].max()))
|
||||
print(f"[fbed] 出图 {time.time()-t0:.2f} s")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
Reference in New Issue
Block a user