Undermind审计前，Cursor重构
This commit is contained in:
Frank14f
@@ -0,0 +1,174 @@
|
||||
# Cylinder benchmark targets for solver validation
|
||||
|
||||
##### [**Undermind**](https://undermind.ai)
|
||||
|
||||
---
|
||||
|
||||
当前阶段不再追求一次覆盖全部圆柱文献，而是先把工作集压缩到两个最有诊断价值的 case：一个固定圆柱文献硬对标，一个旋转圆柱内部回归。这样做的目的是先回答两个更基础的问题：求解器在当前边界能力下是否可靠，以及 `curved_boundary` 这条代码路径是否自洽。
|
||||
|
||||
结合当前代码能力，最合适的固定圆柱主 benchmark 是 \[Sah04\] 的 2D 受限通道圆柱。该 family 与你现有的 `parabolic` 入口和上下 no-slip wall 最一致，且 \[Sah04\] 明确给出了 \\\beta=0.3\\ 下的临界 Reynolds 数和 \\Re=100\\ 时的 Strouhal 数。\[Sah04\] 旋转圆柱则先不做开放来流文献硬对标，因为南北 free-stream 与更匹配的外边界 family 尚未实现；当前最合理的做法，是沿用同一受限通道几何做内部旋转回归，并把 \[Kan99b\] 作为下一阶段开放来流 benchmark 的目标参考。\[Kan99b\]
|
||||
|
||||
在 \[Sah04\] 的表格中，\\\beta=0.30\\ 的第一临界 Reynolds 数和临界 Strouhal 已经给出；下面展示的表格就是当前最适合锁定的文献锚点。上图中的拖曳曲线则说明了为什么 \\\beta=0.30, Re=100\\ 是一个好用的周期态检查点。\[Sah04\]
|
||||
|
||||
## 当前锁定的两个 case
|
||||
|
||||
| ID | 角色 | family | 是否做文献硬评分 | 主要回答的问题 |
|
||||
|:---|:---|:---|:---|:---|
|
||||
| A1 | 固定圆柱主 benchmark | \[Sah04\] confined cylinder | 是 | 现在的 solver 在当前边界能力下能否对标稳定周期 shedding |
|
||||
| B1 | 旋转圆柱内部回归 | internal rotating confined cylinder | 否 | `curved_boundary` 中 moving wall 项，MEA 力，torque，符号约定是否自洽 |
|
||||
|
||||
## A1 固定圆柱主 benchmark
|
||||
|
||||
A1 固定采用 \[Sah04\] 的受限通道 family，并且只选一个最实用的工作点：
|
||||
|
||||
- 阻塞比 \\\beta = D/H = 0.30\\
|
||||
- Reynolds 数 \\Re = 100\\
|
||||
- 入口为 fully developed parabolic inflow \\u=(1-x_2^2,0)\\
|
||||
- 上下边界为 no-slip wall
|
||||
- 圆柱表面为 no-slip curved wall
|
||||
- 出口先用当前最稳的 `neq_extrap`，再用 `zero_gradient` 做敏感性对比
|
||||
- 运行模式统一为 `double_buffer + LES off`
|
||||
|
||||
选择这个点有四个原因。
|
||||
|
||||
- 它和你当前代码能力最匹配 \[Sah04\]
|
||||
- 它已经越过一阶 Hopf 临界点，能同时检查 \\C_d\\、\\C_l\\ 和 \\St\\
|
||||
- \\\beta=0.30\\ 的壁效应足够明显，能放大曲壁与受力统计误差
|
||||
- 文献中给出了明确的 \\Re\_{crit}\\ 和 \\St\\，比只靠图上读趋势更稳 \[Sah04\]
|
||||
|
||||
### A1 的文献定义
|
||||
|
||||
| 项目 | 设定 | 说明 |
|
||||
|:---|:---|:---|
|
||||
| family | `sah04_confined` | 固定圆柱受限通道 |
|
||||
| Reynolds 数定义 | \\Re = U\_{max} D / \nu\\ | 用最大入口速度归一化 \[Sah04\] |
|
||||
| 阻塞比定义 | \\\beta = D/H\\ | \\H\\ 是通道高度 \[Sah04\] |
|
||||
| 入口 | parabolic | \\u=(1-x_2^2,0)\\ \[Sah04\] |
|
||||
| 上下边界 | no-slip wall | 与当前代码能力一致 |
|
||||
| 圆柱边界 | no-slip curved wall | 当前主排错对象 |
|
||||
| 出口 | 先 `neq_extrap` | 文献是二阶导数型出口，你的实现只能近似 |
|
||||
| 上游长度 | 40D | \[Sah04\] |
|
||||
| 下游长度 | 40D | \[Sah04\] |
|
||||
|
||||
### A1 的具体网格方案
|
||||
|
||||
你提到圆柱直径至少要 20 个格点，这个判断是对的。对 \\\beta=0.30\\ 来说，若希望阻塞比在格点上精确，最方便的直径应取 3 的倍数。这样 \\H=D/\beta\\ 才是整数。基于这个原则，当前推荐如下。
|
||||
|
||||
| 方案 | 直径 D | 半径 r | 流体高度 H | 当前代码中的 NY | 流体长度 Lx | 当前代码中的 NX | 圆心 |
|
||||
|:---|---:|---:|---:|---:|---:|---:|:---|
|
||||
| A1 base | 24 | 12.0 | 80 | 82 | 1920 | 1922 | \\(960.5, 40.5)\\ |
|
||||
| A1 confirm | 30 | 15.0 | 100 | 102 | 2400 | 2402 | \\(1200.5, 50.5)\\ |
|
||||
|
||||
这里的换算采用你当前代码的 wall 约定：
|
||||
|
||||
- 流体带高度取 `NY - 2`
|
||||
- 上下边界各占一层边界行
|
||||
- 因此 `NY = H + 2`
|
||||
- 同理，若左右边界节点各占一列，则 `NX = Lx + 2`
|
||||
|
||||
推荐先用 `D=24` 做主排错，因为它已经超过 20 格点，同时保持 \\\beta=0.30\\ 精确。若 A1 结果基本合理，再用 `D=30` 做一次确认，判断误差是边界主导还是分辨率主导。
|
||||
|
||||
### A1 的对标物理量
|
||||
|
||||
A1 不要求所有量都同等对待。当前建议分成三级。
|
||||
|
||||
| 优先级 | 物理量 | 对标方式 | 备注 |
|
||||
|:---|:---|:---|:---|
|
||||
| P1 | \\St\\ | 点目标 | \[Sah04\] 明确给出 \\St=0.2115\\ for \\\beta=0.30, Re=100\\ |
|
||||
| P1 | \\Re\_{crit}\\ | 点目标 | \[Sah04\] 给出 \\Re\_{crit}\approx 94.4\\ for \\\beta=0.30\\ |
|
||||
| P2 | 平均阻力 \\\overline{C_d}\\ | 带宽目标 | 图上可读约 1.8 到 2.0，当前先作为 band target |
|
||||
| P3 | 升力振幅或 RMS | 内部比较 | 文献未在该 case 明确列表，当前先记录，不做硬阈值 |
|
||||
| P3 | 质量漂移 | 诊断量 | 当前不纳入文献分数，但必须长期记录 |
|
||||
|
||||
这意味着：A1 现在最硬的文献锚点其实是 \\Re\_{crit}\\ 和 \\St\\，不是 \\C_l\\ 振幅。平均阻力可以作为第二层目标，但更适合作为带宽检查，而不是一个小数点后三位的硬分数。
|
||||
|
||||
### A1 的通过标准
|
||||
|
||||
| 量 | 通过标准 | 解释 |
|
||||
|:---|:---|:---|
|
||||
| \\St\\ | 相对误差 3 percent 内 | 这是当前最重要的周期态 benchmark |
|
||||
| \\Re\_{crit}\\ | 相对误差 5 percent 内 | 用于确认整体 shedding 触发位置 |
|
||||
| \\\overline{C_d}\\ | 落在 1.8 到 2.0 带内 | 当前只做 band target |
|
||||
| 质量漂移 | 单位 shedding 周期内接近零趋势 | 不要求绝对零，但不能持续单调失控 |
|
||||
|
||||
## B1 旋转圆柱内部回归
|
||||
|
||||
B1 不拿来和 \[Kan99b\] 做硬评分，因为 \[Kan99b\] 用的是开放来流 O-grid 外场与 convective outlet，\[Kan99b\] 而你当前还没有相同 family。B1 的意义不是“验证旋转圆柱文献是否复现”，而是“验证 moving wall 代码路径是否自洽”。
|
||||
|
||||
B1 直接继承 A1 的几何和网格，只改变圆柱壁面速度。这样做能最大限度隔离变量。
|
||||
|
||||
### B1 的统一定义
|
||||
|
||||
| 项目 | 设定 |
|
||||
|:---------------|:------------------------------------|
|
||||
| family | `internal_rotating_confined` |
|
||||
| 几何 | 与 A1 完全相同 |
|
||||
| 参考速度 | 暂统一用 \\U\_{max}\\ |
|
||||
| 自旋比 | \\\alpha = \Omega D / (2U\_{max})\\ |
|
||||
| 入口和上下边界 | 与 A1 完全相同 |
|
||||
| 文献硬评分 | 不做 |
|
||||
|
||||
### B1 的具体工作点
|
||||
|
||||
| 子工况 | 自旋比 \\\alpha\\ | 用途 |
|
||||
|:-------|------------------:|:------------------------------------|
|
||||
| B1-0 | 0.0 | 必须退化回 A1 固定圆柱 |
|
||||
| B1+ | +0.5 | 检查正转时的 moving wall 与受力响应 |
|
||||
| B1- | -0.5 | 检查反转时 lift 和 torque 是否反号 |
|
||||
|
||||
只要 B1+ 和 B1- 不是镜像关系，就先不要扩大到更高 \\\alpha\\ 扫描。
|
||||
|
||||
### B1 的目标物理量
|
||||
|
||||
| 优先级 | 量 | 目标 |
|
||||
|:---|:---|:---|
|
||||
| P1 | \\\overline{C_l}\\ 符号 | 正反转必须反号 |
|
||||
| P1 | torque 符号 | 正反转必须反号 |
|
||||
| P1 | B1-0 与 A1 的一致性 | \\\alpha=0\\ 必须退化成固定圆柱结果 |
|
||||
| P2 | \\\overline{C_d}\\ 随 \\\alpha\\ 的变化 | 作为趋势检查 |
|
||||
| P2 | 质量漂移变化 | 看 moving wall 是否显著放大 leakage |
|
||||
|
||||
## 是否纳入 MRT
|
||||
|
||||
你提到 MRT 可以直接纳入，我同意。MRT 不会改变 benchmark family，本身不会让设置更混乱。当前更合理的做法不是把 MRT 排除，而是把它放在同一 case 的第三个运行层级中。
|
||||
|
||||
| 顺序 | collision | 角色 |
|
||||
|:-----|:----------|:------------------------------------|
|
||||
| 1 | SRT | 最简单可解释基线 |
|
||||
| 2 | TRT | 最直接的边界敏感性对比 |
|
||||
| 3 | MRT | 观察多松弛是否改变 force 和 leakage |
|
||||
|
||||
但解释顺序仍应保持不变：若 SRT 已经明显错，先不要拿 MRT 的较好结果掩盖基础边界问题。
|
||||
|
||||
## 推荐的首轮运行序列
|
||||
|
||||
| 顺序 | run | 目的 |
|
||||
|:-----|:--------------------------------|:--------------------------------------|
|
||||
| 1 | A1 base + SRT + `neq_extrap` | 建立固定圆柱主 benchmark |
|
||||
| 2 | A1 base + TRT + `neq_extrap` | 检查曲壁与碰撞敏感性 |
|
||||
| 3 | A1 base + MRT + `neq_extrap` | 记录 collision family 差异 |
|
||||
| 4 | A1 base + SRT + `zero_gradient` | 判断 outlet 对 \\St\\ 与 force 的影响 |
|
||||
| 5 | A1 confirm + TRT + `neq_extrap` | 做一次分辨率确认 |
|
||||
| 6 | B1-0 + TRT | 检查是否退化回 A1 |
|
||||
| 7 | B1+ + TRT | 检查正转 |
|
||||
| 8 | B1- + TRT | 检查反转 |
|
||||
| 9 | B1+ + MRT | 看 moving wall 下 MRT 是否改变趋势 |
|
||||
|
||||
## 当前阶段的判读重点
|
||||
|
||||
- 若 A1 的 \\St\\ 比 \\\overline{C_d}\\ 更先失真，先查 outlet 和曲壁时序。
|
||||
- 若 A1 的 \\\overline{C_d}\\ 偏差更明显，先查 MEA 力定义、系数归一化和 `q` 分布。
|
||||
- 若 B1 正反转不反号，优先查 `6 w_i (c_i\cdot u_w)` 的符号和 torque 杠杆臂。
|
||||
- 若 B1-0 不能退化回 A1，说明 rotating 路径在 \\\alpha=0\\ 时仍残留额外改动。
|
||||
|
||||
## 结论
|
||||
|
||||
当前最合理的工作集就是 A1 和 B1。A1 负责文献硬对标，且具体锁定在 \[Sah04\] 的 \\\beta=0.30, Re=100\\；网格优先用 `D=24`，再用 `D=30` 做确认。\[Sah04\] B1 负责 rotating path 的内部回归，优先看正反转镜像关系与 torque 符号，而不是先追开放来流文献值。\[Kan99b\] 这样可以在最少的 case 上同时推进“求解器可靠性”和“代码缺陷排查”两条线。
|
||||
|
||||
---
|
||||
|
||||
## References
|
||||
|
||||
\[Sah04\] M. Sahin and R. G. Owens, “A numerical investigation of wall effects up to high blockage ratios on two-dimensional flow past a confined circular cylinder,” Apr. 02, 2004. doi: [10.1063/1.1668285](https://doi.org/10.1063/1.1668285).
|
||||
|
||||
\[Kan99b\] S. Kang, H. Choi, and S. Lee, “Laminar flow past a rotating circular cylinder,” Oct. 07, 1999. doi: [10.1063/1.870190](https://doi.org/10.1063/1.870190).
|
||||
@@ -0,0 +1,258 @@
|
||||
# Boundary features still needed for cylinder validation
|
||||
|
||||
##### [**Undermind**](https://undermind.ai)
|
||||
|
||||
---
|
||||
|
||||
当前代码已经具备做 2D 通道类圆柱验证的主体能力：D2Q9 与 D3Q19，SRT 与 TRT 与 MRT，`double_buffer` 与 `esopull`，可开关 LES，`uniform` 与 `parabolic` 入口，和若干开放出口模式。真正阻塞下一阶段 benchmark 的，不是碰撞模型数量不够，而是少数几个边界功能还没有补齐。它们分别是：南北 free-stream 边界，time-dependent 入口幅值，开放来流 family 的更匹配外边界，以及 z 方向 periodic 边界。\[Qu13, Kan99b, Jia21\]
|
||||
|
||||
这份说明只回答一个问题：这些功能应如何实现，物理上在做什么，数值上应放在什么位置，以及如何兼容你当前的 SRT 与 TRT 与 MRT、`double_buffer` 与 `esopull`、LES 开关。
|
||||
|
||||
## 当前最值得补的功能顺序
|
||||
|
||||
| 顺序 | 功能 | 直接解锁的 benchmark | 为什么先做 |
|
||||
|:---|:---|:---|:---|
|
||||
| 1 | 南北 free-stream 边界 | \[Qu13\] 2D 开放来流 | 数学最简单，直接解锁 2D 恒定来流主 benchmark |
|
||||
| 2 | time-dependent inlet amplitude | \[Joh04\] 全定义复现 | 实现代价低，可完善实现检查 |
|
||||
| 3 | 开放来流 family 的远场与 outlet | \[Kan99b\] 2D 旋转开放来流 | 旋转 benchmark 的关键 |
|
||||
| 4 | z periodic 边界 | \[Jia21\] 3D 固定圆柱 | 一旦实现即可解锁 3D 主 benchmark |
|
||||
|
||||
如果只能先做一个功能，应优先做南北 free-stream 边界。因为它不仅解锁 \[Qu13\]，也会让开放来流中的静止圆柱和旋转圆柱不再被 fixed wall 人为污染。\[Qu13, Kan99b\]
|
||||
|
||||
## 南北 free-stream 边界
|
||||
|
||||
### 物理含义
|
||||
|
||||
\[Qu13\] 的上下边界不是物理壁面，而是远场近似。其假设是：圆柱诱导扰动在足够远的横向距离上已经衰减，因此边界上的流动回到自由来流状态。\[Qu13\]
|
||||
|
||||
对 2D 开放来流圆柱，这意味着在上边界和下边界施加目标宏观状态
|
||||
|
||||
``` math
|
||||
\mathbf{u}_b = (U_{\infty}, 0)
|
||||
```
|
||||
|
||||
``` math
|
||||
\rho_b = \rho_0
|
||||
```
|
||||
|
||||
最实用的实现不是发明一套新的边界家族，而是把你已有的 aligned velocity boundary 泛化到 y 方向边界。
|
||||
|
||||
### 推荐数值形式
|
||||
|
||||
最稳妥的第一版做法是 non-equilibrium extrapolation \[Guo02\]。它的写法是
|
||||
|
||||
``` math
|
||||
f_i(x_b,t+\Delta t)=f_i^{eq}(\rho_b,\mathbf{u}_b)+\left[f_i(x_f,t+\Delta t)-f_i^{eq}(\rho_f,\mathbf{u}_f)\right]
|
||||
```
|
||||
|
||||
这里
|
||||
|
||||
- $`x_b`$ 是边界节点
|
||||
- $`x_f`$ 是相邻内侧流体 donor 节点
|
||||
- 平衡态由目标自由流状态给出
|
||||
- 非平衡部分由 donor 节点外推
|
||||
|
||||
这个形式有几个好处。
|
||||
|
||||
- 与你现有入口边界在代码结构上最接近
|
||||
- 对 SRT 与 TRT 与 MRT 都能共用同一宏观目标状态
|
||||
- 比单纯代数型 local closure 更稳，特别适合工程首版 \[Guo02, Lat08\]
|
||||
|
||||
### 在当前模式中的兼容方式
|
||||
|
||||
| 模块 | 兼容原则 |
|
||||
|:---|:---|
|
||||
| SRT | 直接用当前平衡态与 donor nonequilibrium |
|
||||
| TRT | 同样重构分布，必要时保留现有 `trt_neq_damp` 作为 donor damping |
|
||||
| MRT | 边界先在分布空间重构，再交给常规 MRT bulk 更新 |
|
||||
| LES | 层流 benchmark 默认关。若开启，边界仍使用相同宏观目标，LES 只影响 bulk 局部黏性 |
|
||||
|
||||
### 在 `double_buffer` 与 `esopull` 中的放置位置
|
||||
|
||||
不论 streaming 路径如何，逻辑上都应保持一致：先得到边界节点的 post-stream 已知分布，再重构未知分布。区别只在索引来源。
|
||||
|
||||
| streaming | 实现原则 |
|
||||
|:---|:---|
|
||||
| `double_buffer` | 在 pull streaming 之后，对边界节点用 donor 节点和目标状态重构未知分布 |
|
||||
| `esopull` | 在得到本地 post-stream 视图后，用同一重构公式回填未知方向 |
|
||||
|
||||
因此最好把边界实现写成“给定当前节点、方向集合、目标宏观状态、donor 节点”的统一接口，而不要把公式写死在某一种 streaming 路径中。
|
||||
|
||||
### 工程建议
|
||||
|
||||
第一版不必追求 characteristic far-field。只需先实现 north 与 south 的统一 velocity boundary，支持 `uniform` 目标状态，并把 `boundary_type` 与 `velocity_profile` 解耦。这样同一套 y 边界内核就能同时服务：
|
||||
|
||||
- 开放来流的 free-stream 侧边界
|
||||
- 通道中的移动壁之外的其他速度边界
|
||||
|
||||
## time-dependent inlet amplitude
|
||||
|
||||
### 物理含义
|
||||
|
||||
\[Joh04\] 需要的不是新的空间 profile，而是已有 parabolic profile 的时间幅值调制。\[Joh04\] 因此它是一个很低成本但高收益的功能。
|
||||
|
||||
若把入口 profile 写成模板
|
||||
|
||||
``` math
|
||||
\mathbf{u}_{in}(y,t)=A(t)\,\mathbf{u}_{shape}(y)
|
||||
```
|
||||
|
||||
那么当前已支持的两种 profile 都能自然兼容。
|
||||
|
||||
- `uniform` 对应 $`\mathbf{u}_{shape}=(1,0)`$
|
||||
- `parabolic` 对应 $`\mathbf{u}_{shape}(y)`$ 为固定抛物线形状
|
||||
|
||||
\[Joh04\] 的本质就是给 $`A(t)`$ 一个时间函数。\[Joh04\]
|
||||
|
||||
### 推荐实现
|
||||
|
||||
最简单的工程形式是给 inlet boundary 增加一个可选回调或 schedule。
|
||||
|
||||
| 字段 | 含义 |
|
||||
|:---------------------|:---------------------------|
|
||||
| `profile` | `uniform` 或 `parabolic` |
|
||||
| `amplitude_mode` | `constant` 或 `scheduled` |
|
||||
| `amplitude_value` | 常数幅值 |
|
||||
| `amplitude_schedule` | 按时间步返回幅值的函数或表 |
|
||||
|
||||
这样不需要新边界核，只需要在每一步把目标边界速度更新后交给现有入口重构。
|
||||
|
||||
### 与模式的关系
|
||||
|
||||
它与 SRT 与 TRT 与 MRT、`double_buffer` 与 `esopull`、LES 都无直接耦合。因为变的是目标宏观状态，不是边界算法本身。
|
||||
|
||||
## 开放来流 family 的远场与 outlet
|
||||
|
||||
### 物理含义
|
||||
|
||||
\[Qu13\] 的出口使用 convective boundary。
|
||||
|
||||
``` math
|
||||
\frac{\partial u_i}{\partial t}+C\frac{\partial u_i}{\partial x}=0
|
||||
```
|
||||
|
||||
并取 $`C=1`$。\[Qu13\]
|
||||
|
||||
\[Kan99b\] 也采用了开放外场，并把外边界分成 inflow 半边和 outflow 半边。\[Kan99b\] 对旋转圆柱而言，只修正圆柱壁面速度是不够的；若外场 family 仍不对，升阻力和频率会被远场反射和壁效应同时污染。
|
||||
|
||||
### 对当前工程的建议
|
||||
|
||||
当前已有 `neq_extrap`、`zero_gradient` 和 `blended` 三种 outlet。对于通道类验证，这已经足够先用。但若要更严格地逼近 \[Qu13\] 或 \[Kan99b\]，下一步应考虑增加更接近 convective outlet 的 family。
|
||||
|
||||
最简单的第一版不是 characteristic outlet，而是直接实现一个宏观量层面的 convective update。之后若高 Re 稳定性成为问题，再考虑 characteristic 或 regularized characteristic family \[Izq08, Wis17\]。
|
||||
|
||||
### 推荐路线
|
||||
|
||||
| 阶段 | 功能 |
|
||||
|:-------|:---------------------------------------------|
|
||||
| 第一版 | north 与 south free-stream velocity boundary |
|
||||
| 第二版 | x 出口增加 convective family |
|
||||
| 第三版 | 若高 Re 仍敏感，再加 characteristic outlet |
|
||||
|
||||
### 与模式的关系
|
||||
|
||||
出口 family 应尽量与 collision model 解耦。也就是说：
|
||||
|
||||
- SRT 与 TRT 与 MRT 共享同一个宏观 outlet 目标
|
||||
- 仅 donor nonequilibrium 的处理上保留与现有 config 一致的 damping 参数
|
||||
- LES 不应改变 outlet 类型，只改变内域的有效黏性与非平衡强度
|
||||
|
||||
## z 方向 periodic 边界
|
||||
|
||||
### 物理含义
|
||||
|
||||
\[Jia21\] 的 3D 圆柱 benchmark 之所以能被稳定比较，是因为 spanwise 方向采用 periodic，从而模拟无限长圆柱的一段重复单元。\[Jia21\] periodic 不代表“边界上再做一次边界重构”，而代表拓扑上把两个 z 面缝合成一个连续方向。
|
||||
|
||||
其物理条件是
|
||||
|
||||
``` math
|
||||
\phi(x,y,0,t)=\phi(x,y,L_z,t)
|
||||
```
|
||||
|
||||
对分布函数而言，就是所有跨越 z 边界的 streaming 都直接 wrap 到另一侧。
|
||||
|
||||
### 实现原则
|
||||
|
||||
periodic 边界不属于 bounce-back，也不属于 velocity outlet。它本质上是索引映射。
|
||||
|
||||
对任何具有 $`c_{i,z}=+1`$ 的离散方向，若 pull source 超出上边界，就从另一端取值；反之亦然。
|
||||
|
||||
### 在两种 streaming 路径中的处理
|
||||
|
||||
| streaming | 实现原则 |
|
||||
|:---|:---|
|
||||
| `double_buffer` | pull source 的 z 索引越界时做 modulo wrap |
|
||||
| `esopull` | 针对所有 $`c_{i,z}\neq 0`$ 的方向，读取或写回时做相同 wrap |
|
||||
|
||||
这说明 periodic 的最稳实现位置不在边界核，而在 streaming 索引层。只要 streaming 层完成 wrap，collision、curved boundary、sensor、LES 都不需要知道 periodic 的存在。
|
||||
|
||||
### 与模式的关系
|
||||
|
||||
| 模块 | 关系 |
|
||||
|:------------------|:-------------------------------------------|
|
||||
| SRT 与 TRT 与 MRT | 完全无关，collision 不变 |
|
||||
| LES | 完全无关，局部 SGS 仍按正常 stencil 取邻域 |
|
||||
| 曲壁边界 | 完全无关，圆柱表面仍按现有 kernel 处理 |
|
||||
|
||||
因此 z periodic 是一个高价值、低物理风险的功能。它的难点主要是索引实现，而不是物理公式。
|
||||
|
||||
## 对当前代码架构的建议
|
||||
|
||||
### 先按 boundary family 分层，不要按 case 临时加分支
|
||||
|
||||
建议把缺失功能拆成三层接口。
|
||||
|
||||
| 层 | 作用 |
|
||||
|:---|:---|
|
||||
| `macro target` 层 | 给出边界目标状态，如 $`\rho_b`$ 与 $`\mathbf{u}_b`$ |
|
||||
| `reconstruction` 层 | 用 Zou-He 或 NEQ extrapolation 等方法重构未知分布 |
|
||||
| `streaming topology` 层 | 处理 periodic wrap 或 pull source 定位 |
|
||||
|
||||
这样做的好处是：
|
||||
|
||||
- north 与 south free-stream 只改 `macro target` 与 `reconstruction`
|
||||
- time-dependent inlet 只改 `macro target`
|
||||
- z periodic 只改 `streaming topology`
|
||||
|
||||
### 模式扩展顺序
|
||||
|
||||
建议所有新增功能都按同一顺序落地。
|
||||
|
||||
1. 先在 `double_buffer` 上实现
|
||||
2. 先支持 SRT 基线
|
||||
3. 再验证 TRT 和 MRT
|
||||
4. 最后再移植到 `esopull`
|
||||
5. 层流 benchmark 上默认 LES off
|
||||
|
||||
这是最省调试成本的路线，因为 `double_buffer + SRT + LES off` 最容易隔离边界错误。
|
||||
|
||||
## 最后建议
|
||||
|
||||
如果你的目标是尽快把圆柱验证从当前阶段推进到真正可对标的开放来流和 3D family，那么最值得优先补的不是新的碰撞模型，也不是新的曲壁公式，而是：
|
||||
|
||||
- north 与 south free-stream boundary
|
||||
- inlet amplitude schedule
|
||||
- convective 或更匹配的远场 outlet family
|
||||
- z periodic
|
||||
|
||||
其中第一项和第四项收益最高。第一项直接解锁 \[Qu13\]，第四项直接解锁 \[Jia21\]。\[Qu13, Jia21\] 对旋转圆柱而言，只有在第一项和第三项具备之后，\[Kan99b\] 才能从 internal regression 升级成真正的 benchmark。\[Kan99b\]
|
||||
|
||||
---
|
||||
|
||||
## References
|
||||
|
||||
\[Qu13\] L. Qu, C. Norberg, L. Davidson, S.-H. Peng, and F. Wang, “Quantitative numerical analysis of flow past a circular cylinder at Reynolds number between 50 and 200,” May 01, 2013. doi: [10.1016/J.JFLUIDSTRUCTS.2013.02.007](https://doi.org/10.1016/J.JFLUIDSTRUCTS.2013.02.007).
|
||||
|
||||
\[Kan99b\] S. Kang, H. Choi, and S. Lee, “Laminar flow past a rotating circular cylinder,” Oct. 07, 1999. doi: [10.1063/1.870190](https://doi.org/10.1063/1.870190).
|
||||
|
||||
\[Jia21\] H. Jiang and L. Cheng, “Large-eddy simulation of flow past a circular cylinder for Reynolds numbers 400 to 3900,” Mar. 19, 2021. doi: [10.1063/5.0041168](https://doi.org/10.1063/5.0041168).
|
||||
|
||||
\[Joh04\] V. John, “Reference values for drag and lift of a two‐dimensional time‐dependent flow around a cylinder,” Mar. 10, 2004. doi: [10.1002/FLD.679](https://doi.org/10.1002/FLD.679).
|
||||
|
||||
\[Guo02\] Z. Guo, C. Zheng, and B. Shi, “Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method,” Apr. 01, 2002. doi: [10.1088/1009-1963/11/4/310](https://doi.org/10.1088/1009-1963/11/4/310).
|
||||
|
||||
\[Lat08\] J. Latt, B. Chopard, O. Malaspinas, M. Deville, and A. Michler, “Straight velocity boundaries in the lattice Boltzmann method.” *Physical review. E, Statistical, nonlinear, and soft matter physics*, vol. 77 5 Pt 2, pp. 056703, May 2008, doi: [10.1103/physreve.77.056703](https://doi.org/10.1103/physreve.77.056703).
|
||||
|
||||
\[Izq08\] S. Izquierdo and N. Fueyo, “Characteristic nonreflecting boundary conditions for open boundaries in lattice Boltzmann methods.” *Physical review. E, Statistical, nonlinear, and soft matter physics*, vol. 78 4 Pt 2, pp. 046707, Oct. 2008, doi: [10.1103/PHYSREVE.78.046707](https://doi.org/10.1103/PHYSREVE.78.046707).
|
||||
|
||||
\[Wis17\] G. Wissocq, N. Gourdain, O. Malaspinas, and A. Eyssartier, “Regularized characteristic boundary conditions for the Lattice-Boltzmann methods at high Reynolds number flows,” *J. Comput. Phys.*, vol. 331, pp. 1–18, Jan. 2017, doi: [10.1016/j.jcp.2016.11.037](https://doi.org/10.1016/j.jcp.2016.11.037).
|
||||
File diff suppressed because one or more lines are too long
+628
@@ -0,0 +1,628 @@
|
||||
Physics of Fluids
|
||||
[Non-Text]
|
||||
AIP Publishing
|
||||
09 May 2026 02:44:09
|
||||
RESEARCH ARTICLE | MAY 01 2004
|
||||
# A numerical investigation of wall effects up to high blockage ratios on two-dimensional flow past a confined circular cylinder
|
||||
Mehmet Sahin; Robert G. Owens
|
||||
|
||||
|
||||
Check for updates
|
||||
Physics of Fluids 16, 1305–1320 (2004)
|
||||
https://doi.org/10.1063/1.1668285
|
||||
|
||||
|
||||
|
||||
View Online
|
||||
|
||||
|
||||
|
||||
|
||||
Export Citation
|
||||
|
||||
# Articles You May Be Interested In
|
||||
Shear-induced autorotation of freely rotatable cylinder in a channel flow at moderate Reynolds number
|
||||
Physics of Fluids (April 2018)
|
||||
Frequency lock-in mechanism in the presence of blockage effects
|
||||
Physics of Fluids (July 2024)
|
||||
The blockage and erosion characteristics of woody debris flow on an erodible gully bed: Insight from a small-scale model experiment
|
||||
Physics of Fluids (December 2024)
|
||||
|
||||
|
||||
# AIP Advances
|
||||
# Why Publish With Us?
|
||||
|
||||
|
||||
|
||||
21DAYS average time to1st decision
|
||||
|
||||
|
||||
|
||||
|
||||
OVER 4 MILLION views in the last year
|
||||
|
||||
|
||||
|
||||
|
||||
INCLUSIVE scope
|
||||
|
||||
Learn More
|
||||
|
||||
|
||||
|
||||
AIP Publishing
|
||||
|
||||
09 May 2026 02:44:09
|
||||
HTmL AB:STRACT * LINKS
|
||||
PHYSICS OF FLUIDS
|
||||
VOLUME 16, NUMBER 5
|
||||
MAY 2004
|
||||
# A numerical investigation of wall effects up to high blockage ratios on two-dimensional flow past a confined circular cylinder
|
||||
Mehmet Sahin and Robert G. Owensa)
|
||||
LMF-ISE-FSTI, Ecole Polytechnique Fe´de´rale de Lausanne, CH 1015 Lausanne, Switzerland
|
||||
~Received 25 September 2003; accepted 20 January 2004; published online 2 April 2004!
|
||||
A finite volume method based on a velocity-only formulation is used to solve the flow field around a confined circular cylinder in a channel in order to investigate lateral wall proximity effects on stability, Strouhal number, hydrodynamic forces and wake structure behind the cylinder for a wide range of blockage ratios (0.1,b<0.9) and Reynolds numbers $( 0 < R e \leqslant 2 8 0 )$ . For blockage ratios less than approximately 0.85 a first critical Reynolds number is identified at which a supercritical Hopf bifurcation of the symmetric solution occurs. For blockage ratios greater than about 0.687 and at Reynolds numbers exceeding the first critical Reynolds number a second curve of neutral stability is seen, representing a pitchfork bifurcation of the steady symmetric solution to one of two possible steady asymmetric solutions. Either side of the neutral stability curve for the pitchfork bifurcation our linear stability analysis and direct numerical simulations demonstrate that although the flow is linearly stable it is unstable to finite two-dimensional perturbations. At blockage ratios larger than about 0.82 the steady asymmetric solutions also become unstable through a Hopf bifurcation. In contrast with the first Hopf bifurcation of the symmetric solution at lower Reynolds numbers numerical calculations of the lift coefficient reveal that the oscillations are no longer symmetric in the rising and falling parts of each cycle. Very strong vortices shed from the cylinder and the wall cause drastic increases in the amplitudes of the lift and drag coefficients. A co-dimension 2 point where pitchfork and Hopf bifurcations occur simultaneously has been located in parameter space. Altogether, four distinct regions in the parameter space $( \beta , R e ) \in ( 0 , 0 . 9 ] \times ( 0 , 2 8 0 ]$ have been identified, each corresponding to a different class of flow: ~i! Steady symmetric flow, ~ii! symmetric vortex shedding, ~iii! steady asymmetric flow, and ~iv! asymmetric vortex shedding, where a periodic-in-time flow is classed as symmetric or asymmetric depending on whether the time-average over one cycle of the lift coefficient is zero or not. Numerical solutions are computed on meshes having up to 1.8 million degrees of freedom. Extensive comparisons are made with the results available in the literature. © 2004 American Institute of Physics. @DOI: 10.1063/1.1668285#
|
||||
# I. INTRODUCTION
|
||||
It is no exaggeration to say that an enormous ~and still rapidly growing! corpus of literature on the subject of bluff body wakes has developed since the pioneering work of von Ka´rma´n early last century. This fact is an attestation to both the difficulty in understanding and adequately describing the flow bifurcations that occur at various values of the Reynolds number in viscous flows and the interest in doing so. Flows having particularly simple setups such as those past a sphere or cylinder have succeeded in drawing experimentalists, theoreticians and computational fluid dynamicists into the fray that has gone on through the decades and only very recently are consensuses emerging.
|
||||
Details of recent theoretical, experimental and computational developments for unbounded flow past a cylinder may be found in the review paper of Williamson,1 where particular attention is paid to the vortex dynamics in the cylinder wake. Our interest in this paper is a careful analysis of lateral wall effects on viscous flow past a confined cylinder. What we have in mind is depicted in Fig. 1. In this figure an infinitely long cylinder of diameter D is placed symmetrically between parallel lateral walls a distance H apart. The parameter $\beta { \equiv } D / H$ is usually termed the blockage ratio. In stark contrast to the wealth of insight and commentary available on vortex dynamics in the wake of an unbounded cylinder we find ourselves with only a handful of papers offering a serious treatment of the blockage ratio effects present in the confined cylinder problem. This paucity of scientific literature should not be interpreted as implying that the problem is an unimportant one, however. On the contrary, even for unbounded flow past a cylinder the ~infinite! flow domain has to be replaced with ~or mapped onto! a finite one, thus introducing numerical or experimental blockage effects that may have considerable influence over the determined values of the flow parameters.2,3 Many of the blockage ratio effects described in the literature are more or less evident:
|
||||
|
||||
~1! In the steady flow regime, bringing the walls closer to the cylinder results in the appearance of the twin vortices in the cylinder wake at higher Reynolds numbers.
|
||||
~2! At any given modest (&50) Reynolds number and for
|
||||
a! Author to whom correspondence should be addressed. Electronic mail: robert.owens@epfl.ch
|
||||
1070-6631/2004/16(5)/1305/16/$22.00
|
||||
1305
|
||||
© 2004 American Institute of Physics
|
||||
09 May 2026 02:44:09
|
||||
Phys. Fluids, Vol. 16, No. 5, May 2004
|
||||
M. Sahin and R. G. Owens
|
||||
|
||||
|
||||
|
||||
FIG. 1. Schematic of a cylinder placed symmetrically in a plane channel. The cylinder diameter is D and the channel height H.
|
||||
|
||||
$\beta { \leqslant } 0 . 2$ the length of the closed vortex bubble decreases with wall proximity, while remaining a linear function of R e . 2– 4 $R e . ^ { 2 - 4 }$
|
||||
~3! For increasing blockage ratios $\beta$ up to 0.5 the steady two-dimensional base flow is stabilized with respect to infinitesimal perturbations due to constraint by the confining walls of the separating shear layer that exists between the cylinder wake and the wall boundary layer vorticity.3,4
|
||||
~4! Once the critical Reynolds number for the primary instability has been exceeded the frequency with which periodic two-dimensional vortex shedding takes place at a given Reynolds number is an increasing function of $\bar { \boldsymbol { \beta } } . ^ { 2 , 4 - 6 }$ ~Note that the spurious result obtained by Stansby and Slaouti6 for $\beta = 0 . 5$ is thought to be due to neglect of the boundary layers in their numerical simulation using random vortex methods.!
|
||||
~5! At $R e = O ( 1 0 0 )$ both the mean drag coefficient $C _ { d }$ and the separation angle of the vortex bubble increase as the walls approach the cylinder.2,3,5,6
|
||||
In addition to the obvious interest of wall blockage effects and as observed by Chen $e t a l . , { ^ 4 }$ the choice of a bounded domain allows a more definitive specification of the flow ~both numerically and experimentally! than is possible in the unbounded case, whilst conserving the essential features of the latter.
|
||||
A problem bearing some similarities to that of flow past a confined cylinder is that of flow around a cylinder placed at various heights above a plane boundary. A recent literature survey of experimental investigations into this problem may be found in the paper of Lei et al.7 These studies have sought to address the issue of how forces on the cylinder and vortex shedding frequency depend on the ratio $g / D$ of the gap between the cylinder and the wall, g, and the cylinder diameter, D. They have also been concerned with understanding the effect on these quantities of the boundary layer thickness and the velocity gradient. Most of the experiments have been conducted at Reynolds numbers in the sub-critical regime $\left[ R e = O ( 1 \times 1 0 ^ { 4 } ) \right]$ in which the boundary layer is still laminar. Lei et $a l . ^ { 7 }$ found that the drag coefficient $C _ { d }$ increased with increasing gap ratio because of the reduction in the base pressure. The same trend in base pressure dependence had been observed by Bearman and Zdravkovich.8 The latter authors further found that the Strouhal number for $g / D { \gtrsim } 0 . 3$ was more or less constant in their experiments at a Reynolds number of $4 . 8 \times 1 0 ^ { 4 }$ . Lei et $a l . ^ { 7 }$ also noted only slight fluctuations in a Strouhal number computed from the free-stream velocity for a similar range of gap to diameter ratios. However, for gaps less than 0.3 cylinder diameters8 or $0 . 2 \mathrm { - } 0 . 3$ diameters ~depending on the boundary layer thickness7 !, vortex shedding was suppressed. Differences in the quantification of the vortex shedding suppression gap ratio were due possibly to differences in the boundary layer thicknesses generated by the experimentalists and also to the manner in which the critical gap ratio was identified: Bearman and $\mathrm { Z d r a v k o v i c h } ^ { 8 }$ using a spectral analysis of hot-wire signals in the cylinder wake whereas the method of Lei $e t a l . ^ { 7 }$ was based on observation of the spectrum of the lift coefficient. Suppression of vortices for a sufficiently small gap ratio was also confirmed by Zovatto and Pedrizzetti,9 who used a finite element method based on a vorticity-streamfunction formulation to analyze flow around a cylinder positioned eccentrically between two lateral walls. For very small gap ratios Zovatto and Pedrizzetti9 found a recirculating bubble on the wall downstream of the cylinder. A separation bubble on the wall had also been seen earlier by Bearman and Zdravkovich8 for gap ratios smaller than the critical value for vortex suppression. We will return to some of these flow phenomena in our discussion of our numerical results in Sec. IV for large blockage ratios.
|
||||
|
||||
The motivation for the present study is twofold. First, the rich fluid dynamics in the wake and near the lateral walls deserves to be investigated with greater numerical accuracy than has been possible with the computational resources available to other researchers at the time at which they prepared their manuscripts. Computations on meshes allowing for only tens of thousands of degrees of freedom have been typical ~for example, Refs. 4 and 10!. In the present study a novel finite volume method $\begin{array} { r l } { { 1 1 - 1 3 } } & { { } } \\ { . } \end{array}$ is used in a parallel implementation, permitting up to 1.8 million degrees of freedom and thus a higher resolution of the wake and boundary layer structures. Second, the only previous numerical linear stability analysis of flow past a confined cylinder available to $\mathrm { u s } ^ { 4 }$ went no further than a blockage ratio of $\beta { = } 0 . 7$ . From the results of this publication the trend seemed to be one of decreasing linear stability of the two-dimensional flow for $\beta { > } 0 . 5$ . Stability was always lost over the range of blockage ratios considered through a symmetry-breaking supercritical Hopf bifurcation. We wish in this paper to investigate the effect on the critical Reynolds number of choosing $\beta { > } 0 . 7$ and to identify the nature of the flow instabilities by means of an Arnoldi method.
|
||||
The outline of the present paper is as follows: In Sec. II we describe the problem to be solved and furnish the reader with a brief description of the numerical method used to analyze flow past a confined cylinder at Reynolds numbers up to 280. Section III is dedicated to validation of our numerical scheme for the classical problem of unbounded twodimensional flow past a circular cylinder. Extensive comparison with other results in the literature is made. In particular we find excellent agreement with previously obtained values for the drag coefficients, first critical Reynolds numbers $R e _ { \mathrm { c r i t 1 } }$ and the corresponding critical Strouhal numbers. In
|
||||
1306
|
||||
09 May 2026 02:44:09
|
||||
Phys. Fluids, Vol. 16, No. 5, May 2004
|
||||
A numerical investigation of wall effects
|
||||
Sec. IV we are concerned with a detailed description of wake dynamics and interactions of the wake and wall boundary layers for blockage ratios up to 0.9. For blockage ratios below approximately 0.85 the locus of a supercritical Hopf bifurcation may be traced out in parameter space. At higher Reynolds numbers and for blockage ratios sufficiently large there is a pitchfork bifurcation of the steady symmetric state to one of two asymmetric steady states. Either side of the curve of neutral stability for the pitchfork bifurcation the steady solutions are linearly stable but appear on the basis of direct numerical simulations to be unstable to finite twodimensional perturbations. For yet larger Reynolds numbers and $\beta { \gtrsim } 0 . 8 2$ a Hopf bifurcation of the asymmetric state occurs. The oscillations are now quite different from those associated with the first symmetry-breaking instability, the amplitude of the drag and lift coefficients being much stronger and the oscillations are now asymmetric in time. Finally, we draw some conclusions.
|
||||
# II. MATHEMATICAL PROBLEM AND NUMERICAL SCHEME
|
||||
An infinitely long cylinder of diameter D is placed midway between two parallel planes which are a distance H apart, as shown in Fig. 1. Let us denote by $U _ { \mathrm { m a x } }$ the maximum inlet fluid speed. The incompressible unsteady Navier– Stokes equations may be written in dimensionless form as
|
||||
$$
|
||||
\frac {\partial \mathbf {u}}{\partial t} + (\mathbf {u} \cdot \nabla) \mathbf {u} = - \nabla p + \frac {1}{R e} \nabla^ {2} \mathbf {u}, \tag {1}
|
||||
$$
|
||||
$$
|
||||
\nabla \cdot \mathbf {u} = 0, \tag {2}
|
||||
$$
|
||||
where, in the usual notation, $\mathbf { u } { = } ( u _ { 1 } , u _ { 2 } )$ denotes the velocity field, p the pressure and Re is a Reynolds number. In the present work the Reynolds number is defined as Re $= U _ { \mathrm { m a x } } D / v$ where v is the kinematic viscosity. In the presentation of results in Secs. III and IV for those flows exhibiting periodic vortex shedding, a Strouhal number St is defined by $S t { = } D / ( T U _ { \operatorname* { m a x } } ) .$ , where T is the period of vortex shedding. We denote by $( \mathbf { x } , t ) { = } ( ( x _ { 1 } , x _ { 2 } ) , t )$ a generic point in space and time.
|
||||
In Sec. III we approximate the unbounded cylinder geometry by choosing $\beta = 0 . 0 1$ and the following boundary conditions:
|
||||
$\mathrm { C y l i n d e r ~ s u r f a c e : } ~ { \bf u } { = } ( 0 , 0 ) ,$ ~3!
|
||||
$\mathrm { L a t e r a l ~ w a l l s } \colon \mathrm { ~ } \mathbf { u } = ( 1 , 0 ) ,$ ~4!
|
||||
$\mathrm { I n f l o w : } \qquad \mathbf { u } = ( 1 , 0 ) ,$ ~5!
|
||||
$\mathrm { O u t f l o w : } \qquad { \frac { \partial ^ { 2 } u _ { 1 } } { \partial x _ { 1 } ^ { 2 } } } = 0 , \ { \frac { \partial u _ { 2 } } { \partial x _ { 1 } } } = 0 .$ 9u1 du2 ~6! x1
|
||||
For the confined cylinder problem ~see Sec. IV! Eqs. ~1! and ~2! are solved subject to the following boundary conditions on the components of velocity:
|
||||
$\mathrm { C y l i n d e r ~ s u r f a c e : } ~ { \bf u } { = } ( 0 , 0 ) ,$ ~7!
|
||||
$\mathrm { L a t e r a l ~ w a l l s } \colon \mathrm { ~ } \mathbf { u } = ( 0 , 0 ) ,$ ~8!
|
||||
$\mathrm { I n f l o w : } \qquad \mathbf { u } = ( 1 - x _ { 2 } ^ { 2 } , 0 ) ,$ ~9!
|
||||
$$
|
||||
\text { Outflow: } \quad \frac {\partial^ {2} u _ {1}}{\partial x _ {1} ^ {2}} = 0, \frac {\partial u _ {2}}{\partial x _ {1}} = 0. \tag {10}
|
||||
$$
|
||||
Some care needs to be taken with how the second normal derivative outflow condition is imposed, due to possible linear dependence in the discrete equation set of this condition with the discrete form of the continuity equation. More precisely stated, the second normal derivative outflow condition will be automatically satisfied at $x _ { 2 } = 0$ upon imposition of the continuity equation ~2! within each finite volume. For all the results presented in Sec. III the dimensionless upstream and downstream channel lengths were set equal to 400 D. In Sec. IV these lengths were chosen to both be equal to 40 D. The choice of outflow boundary conditions ~6! and ~10! was motivated by the fact that our numerical method uses a velocity-only formulation so that the usual traction-free conditions could not easily be implemented. The free boundary layer type of conditions ~6! and ~10! were used successfully by Kourta et al.14 in finite volume simulations of a twodimensional plane mixing layer. Although Jin and Braza15 later developed a nonreflecting outlet condition that greatly reduced feedback noise when compared with the outlet condition of Kourta et al.,14 the outlet length used for the calculations in the present paper are considered sufficiently great that the difference between the influence of the one set of exit conditions and the other on drag, linear stability and Strouhal number would be negligible. The more complicated exit conditions of Jin and Braza15 are therefore not implemented.
|
||||
Let n denote a unit outward pointing normal vector to the boundary ]V of a finite volume V. Then integration of ~2! over V and taking the vector product of ~1! with n, followed by integration around ]V leads, respectively, to
|
||||
$$
|
||||
\oint_ {\partial \Omega} \mathbf {n} \cdot \mathbf {u} d s = 0 \tag {11}
|
||||
$$
|
||||
and
|
||||
$$
|
||||
\oint_ {\partial \Omega} \mathbf {n} \times \left[ \frac {\partial \mathbf {u}}{\partial t} + (\nabla \times \mathbf {u}) \times \mathbf {u} + \frac {1}{R e} \nabla \times (\nabla \times \mathbf {u}) \right] d s = \mathbf {0}. \tag {12}
|
||||
$$
|
||||
In our numerical scheme the continuity equation ~11! is satisfied within each finite volume while ~12! is applied to each finite volume except the finite volumes next to the wall. Therefore, vorticity creation is allowed within these finite volumes in order to satisfy the no-slip boundary conditions. Equations ~11! and ~12! with no-slip boundary conditions are enough to solve the problem in a simply connected domain ~such as that found in the lid-driven cavity problem, for $\mathrm { e x a m p l e } ^ { 1 2 } )$ . However, if the domain is not simply connected there is a need for additional equations. This is because there is a potential problem in our velocity-only formulation with multi-valuedness of the pressure field, even though the pressure does not appear explicitly as a dependent variable in our formulation. To rectify this a Kutta-type condition
|
||||
$$
|
||||
\oint_ {\Gamma} \mathbf {n} \times \left[ \frac {\partial \mathbf {u}}{\partial t} + (\nabla \times \mathbf {u}) \times \mathbf {u} + \frac {1}{R e} \nabla \times (\nabla \times \mathbf {u}) \right] d s = \mathbf {0}, \tag {13}
|
||||
$$
|
||||
1307
|
||||
09 May 2026 02:44:09
|
||||
Phys. Fluids, Vol. 16, No. 5, May 2004
|
||||
M. Sahin and R. G. Owens
|
||||
|
||||
TABLE I. Values of grid parameters $i _ { \mathrm { m a x } } , k _ { \mathrm { m a x } } , k _ { \mathrm { w a l l } }$ , and N .
|
||||
|
||||
<table><tr><td rowspan="2"><eq>\beta</eq></td><td colspan="4">M1</td><td colspan="4">M2</td><td colspan="4">M3</td></tr><tr><td><eq>i_{\text{max}}</eq></td><td><eq>k_{\text{max}}</eq></td><td><eq>k_{\text{wall}}</eq></td><td>N</td><td><eq>i_{\text{max}}</eq></td><td><eq>k_{\text{max}}</eq></td><td><eq>k_{\text{wall}}</eq></td><td>N</td><td><eq>i_{\text{max}}</eq></td><td><eq>k_{\text{max}}</eq></td><td><eq>k_{\text{wall}}</eq></td><td>N</td></tr><tr><td>0.01</td><td>181</td><td>301</td><td>137</td><td>89 336</td><td>361</td><td>601</td><td>273</td><td>355 312</td><td>721</td><td>1201</td><td>545</td><td>1 417 184</td></tr><tr><td>0.1</td><td>181</td><td>441</td><td>77</td><td>116 536</td><td>361</td><td>881</td><td>153</td><td>462 512</td><td>721</td><td>1761</td><td>305</td><td>1 842 784</td></tr><tr><td>0.2</td><td>181</td><td>421</td><td>57</td><td>109 336</td><td>361</td><td>841</td><td>113</td><td>433 712</td><td>721</td><td>1681</td><td>225</td><td>1 727 584</td></tr><tr><td>0.3</td><td>181</td><td>411</td><td>47</td><td>105 736</td><td>361</td><td>821</td><td>93</td><td>419 312</td><td>721</td><td>1641</td><td>185</td><td>1 669 984</td></tr><tr><td>0.5</td><td>181</td><td>401</td><td>37</td><td>102 136</td><td>361</td><td>801</td><td>73</td><td>404 912</td><td>721</td><td>1601</td><td>145</td><td>1 612 384</td></tr><tr><td>0.7</td><td>181</td><td>391</td><td>27</td><td>98 536</td><td>361</td><td>781</td><td>53</td><td>390 512</td><td>721</td><td>1561</td><td>105</td><td>1 554 784</td></tr><tr><td>0.9</td><td>181</td><td>381</td><td>17</td><td>94 936</td><td>361</td><td>761</td><td>33</td><td>376 112</td><td>721</td><td>1521</td><td>65</td><td>1 497 184</td></tr></table>
|
||||
is imposed around the closed path G formed from the union of the outer edges of the finite volumes on the cylinder surface. The condition ~13! guarantees that
|
||||
$$
|
||||
\oint_ {\Gamma} \mathbf {n} \times \nabla p d s = \mathbf {k} [ p ] = \mathbf {0}, \tag {14}
|
||||
$$
|
||||
where k is a unit vector normal to the plane of the flow, and @ p# denotes the jump in the pressure on passing once around G. Since ~12! is satisfied in every interior finite volume, satisfaction of ~13! ensures that $p$ is single-valued at every interior finite volume vertex. The pressure can be obtained by integrating the two components of the pressure gradient appearing in the equations of linear momentum in a manner analogous to that used in finding a streamfunction from a given velocity. The values of p on the domain boundaries, when required, are determined by first computing ]p/]n from ~1!.
|
||||
A fully implicit second-order cell-vertex finite volume method based on a velocity-only formulation is used for the discretization of ~1! and ~2!. Discretization of the integrals appearing in ~11! and ~12! is effected by using the mid-point rule on cell faces. Full details of the method are supplied in two recent papers by the present authors11,12 and, in the interests of brevity, will not be reproduced here. For the timedependent computations presented in Secs. III and IV we discretize in time using an Euler implicit method and for computing steady-state base flows a Newton method is employed.
|
||||
A major part of the present paper is concerned with the linear stability of two-dimensional flow at various different blockage ratios. Consider the perturbed flow
|
||||
$$
|
||||
\mathbf {u} (\mathbf {x}, t) = \mathbf {U} (\mathbf {x}) + \mathbf {v} (\mathbf {x}) \exp (\sigma t), \tag {15}
|
||||
$$
|
||||
where U(x) is the ~numerically determined! steady base flow at a given Reynolds number. Then discretizing the dimensionless Navier–Stokes equations as described above leads to an algebraic system of equations
|
||||
$$
|
||||
\mathbf {A} \mathbf {x} = \sigma \mathbf {M} \mathbf {x}, \tag {16}
|
||||
$$
|
||||
for the nodal values of the perturbation velocity v. The matrices A and M in ~16! are block quad-diagonal and block bi-diagonal, respectively. The GEVP ~16! may be solved by applying Arnoldi’s method16,17 to the equivalent system
|
||||
$$
|
||||
\mathbf {C} \mathbf {x} = \mu \mathbf {x}, \tag {17}
|
||||
$$
|
||||
where $\mathbf { C } = ( \mathbf { A } - \lambda \mathbf { M } ) ^ { - 1 } \mathbf { M }$ and $\mu { = } \left( \sigma { - } \lambda \right) ^ { - 1 }$ . Application of the Arnoldi method results in the construction of an upper Hessenberg matrix whose eigenvalues are approximations to a subset of the eigenvalues $\mu$ of C. From the properties of Arnoldi’s method and in the absence of a shift $\lambda , ^ { 1 8 }$ best resolution of the s-spectrum is expected to be near the origin.
|
||||
The coefficient matrix A in ~16! is almost identical ~by construction! to that which arises in the computations of the steady base flow using Newton’s method. Solutions to all the discrete algebraic equations that arise in the steady, unsteady or eigenvalue problems of this paper have been obtained by implementing the MUltifrontal Massively Parallel Solver ~MUMPS! of Amestoy et al.19,20 The multifrontal method used is a direct method based on LU decomposition for the solution of sparse systems of linear equations with optimum fill in. The algorithms employed by MUMPS use a dynamic distributed task scheduling technique that permits numerical pivoting and the transfer of computational tasks to lightly loaded processors. The calculations have been performed on an SGI Origin 3800 parallel machine with 124 processors and on a Linux cluster with 22 processors.
|
||||
Three different finite volume grids (M 1 – M 3) have been used for each value of the blockage ratio considered in this paper. Each of the meshes has been generated algebraically and then smoothed by solving elliptic partial differential equations for the spatial variables $x _ { 1 }$ and $x _ { 2 }$ where derivatives are with respect to mapped variables in a space in which the mesh appears rectangular.21 For the present problem the physical grid is ‘‘cut’’ along the line $x _ { 2 } = 0$ from the rear stagnation point to the outlet before being mapped. The method of Steger and Sorenson21 allows both grid cell sizes and grid cell skewness to be controlled at the inner and outer boundaries. Meshes M 1 to M 3 are characterized by $i _ { \mathrm { m a x } }$ : The number of nodes on the surface of the cylinder, $k _ { \mathrm { m a x } } \colon$ The number of nodes along the line $x _ { 2 } = 0$ from the rear stagnation point on the cylinder to the outflow boundary and $k _ { \mathrm { w a l l } }$ : The number of nodes in the gap between the cylinder and a lateral wall. The values of $i _ { \mathrm { m a x } } , k _ { \mathrm { m a x } } , k _ { \mathrm { w a l l } }$ , and N ~the number of degrees of freedom! for the three meshes are supplied in Table I for different blockage ratios.
|
||||
# III. FLOW PAST AN UNBOUNDED CIRCULAR CYLINDER
|
||||
Flow around an unbounded circular cylinder is a classical benchmark problem for which a large number of numerical and experimental results exists. In this problem, and in approximation to the case of an unbounded flow domain, a circular cylinder of diameter D51.00 is placed symmetrically in a channel with blockage ratio $\beta = 0 . 0 1$ . For the numerical linear stability analysis the three meshes M 1 to M 3 were used, with $i _ { \mathrm { m a x } } , \ : k _ { \mathrm { m a x } }$ , and $k _ { \mathrm { w a l l } }$ as given in Table I. However, for unsteady time-dependent simulations we were only able to afford to use M1 and M2, the unsteady calculations on M3 proving to be prohibitively expensive. On the boundaries of the computational domain the conditions ~3!– ~6! were imposed.
|
||||
1308
|
||||
09 May 2026 02:44:09
|
||||
Phys. Fluids, Vol. 16, No. 5, May 2004
|
||||
A numerical investigation of wall effects
|
||||
|
||||
TABLE II. Unbounded flow past a cylinder. Comparison of critical Reynolds numbers computed on M1–M3 with others in the literature.
|
||||
|
||||
<table><tr><td></td><td>M1</td><td>M2</td><td>M3</td><td>Extrapolated</td><td>Jackson</td><td>Ding and Kawahara</td><td>Noack and Eckelmann</td><td>Chen et al.</td></tr><tr><td><eq>Re_{\text{crit1}}</eq></td><td>47.08</td><td>46.82</td><td>46.76</td><td>46.74</td><td>46.184</td><td>46.389</td><td>50</td><td>47.9</td></tr><tr><td><eq>St_{\text{crit1}}</eq></td><td>0.1163</td><td>0.1166</td><td>0.1167</td><td>0.1167</td><td>0.13804</td><td>0.12619</td><td>0.132</td><td>0.138</td></tr></table>
|
||||
|
||||
The linear stability analysis predictions of the critical Reynolds and Strouhal numbers corresponding to the onset of the first flow instability are supplied in Table II, as computed on meshes M 1 – M 3. Also shown are the values of these quantities when extrapolated to zero mesh size. The extrapolated critical Reynolds number is found to be $R e _ { \mathrm { c r i t 1 } }$ 546.74 with a corresponding Strouhal number of $S t _ { \mathrm { c r i t 1 } }$ 50.1167. These values are compared with others in the literature in the same table. Although we find good agreement for the critical Reynolds number with the result of Ding and Kawahara22 $( R e _ { \mathrm { c r i t 1 } } = 4 6 . 3 8 9 )$ , and Jackson10 $( R e _ { \mathrm { c r i t 1 } }$ 546.184), the critical Strouhal number manifests wider scatter in the cited references. Issues such as the blockage ratios chosen, distances from the cylinder of the upstream and downstream boundaries, boundary conditions, mesh resolution and number of eigenvalues determined may be amongst the reasons for discrepancies in the numerical results. In addition to our mesh convergence study, we present a convergence study of the leading eigenvalues on mesh M 1 with the Krylov subspace dimension m and shift parameter l in Table III in order to show that our leading eigenvalue is essentially independent of both m and l for sufficiently large values of these two parameters. Although the leading eigenvalue converges very rapidly with a suitably chosen complex shift around the leading eigenvalue, it requires complex arithmetic. A real shift also dramatically improves the convergence of the leading eigenvalue while avoiding complex arithmetic which significantly increases the memory requirements during LU factorization. Our calculations show that a Krylov subspace dimension as low as 250 can be enough to compute the leading eigenvalue with l50.5060.00i while with no shift a Krylov subspace dimension larger than 1000 may be required.
|
||||
|
||||
The computed eigenspectrum on mesh M 3 at the critical Reynolds number is given in Fig. 2. Although we present the first 250 computed eigenvalues, calculations with higher Krylov subspace dimensions showed that only the leading eigenvalues and the eigenvalues around the origin were properly converged. As may be seen, the most dangerous eigenvalue pair is well separated from the rest of the spectrum, unlike the eigenspectrum for the two-dimensional liddriven cavity problem, for example.13 This is likely to be the reason for well-developed periodic flow observed far beyond the critical Reynolds number. Our critical Strouhal number ~0.1167! compares very well with the Strouhal number St 50.1179 computed at the same Reynolds number ~46.74! from a curve fit of the two-dimensional experimental data of Williamson.23 In addition, our critical Strouhal number agrees quite well with the Strouhal number $( S t = 0 . 1 1 8 3 4 )$ of the direct numerical simulation of Posdziech and Grundmann,24 even though their critical Strouhal number and ours were computed at two slightly different Reynolds numbers $( R e = 4 7 . 5 0$ and $R e = 4 6 . 7 4 ,$ , respectively!.
|
||||
In Fig. 3 we present comparisons of the Strouhal number versus Reynolds number and in Fig. 4 comparisons of the drag coefficient $C _ { d } { = } F _ { x } / 0 . 5 U _ { \mathrm { m a x } } ^ { 2 } D$ versus Reynolds number, in further verification of our numerical scheme. Our Strouhal numbers are seen to be in very good agreement with those from the experimental work of Williamson23 for Reynolds numbers up to 200. Beyond this point the flow becomes three-dimensional and we do not expect to have agreement with the experimental results. Good agreement for St and Cd with results from the two-dimensional numerical simulations of Henderson25 and Posdziech and Grundmann24 may also be seen from Figs. 3 and 4. Our computed lift coefficients $C _ { l }$ at Reynolds numbers of 100 and 200 are 60.3333 and 60.6861 and these are in satisfactory agreement with Posdziech and Grundmann’s values of 60.321 04 and 60.673 15, respectively. Although both Henderson25 and Posdziech and Grundmann24 used high-order spectral elements the differences between their two sets of results are due to the use of different blockage ratios in their calculations. However, as the Reynolds number increases the difference in their computed results becomes smaller. An interesting convergence study on the extension of the computational domain boundary is given by Posdziech and Grundmann24 at $R e = 2 0 0 . 0 0$ . The authors concluded that the lateral boundaries should be set at a distance of at least 70 diameters away in order to obtain a Strouhal number independent of yet smaller blockage ratios. At lower Reynolds number the effect of the lateral boundaries becomes more severe ~see Fornberg,26 for example!. In addition, Zisis and Mitsoulis27 showed that the convergence of the total drag at $R e { = } 0 . 0 0$ may be very poor as $\beta$ goes to zero.
|
||||
|
||||
TABLE III. Unbounded flow past a cylinder. Convergence of the leading eigenvalue at Re547.08 on M1 with the Krylov space dimension m and shift parameter l.
|
||||
|
||||
<table><tr><td>m</td><td><eq>\lambda = 0.00 \pm 0.00i</eq></td><td><eq>\lambda = 0.25 \pm 0.00i</eq></td><td><eq>\lambda = 0.50 \pm 0.00i</eq></td></tr><tr><td>250</td><td><eq>-2.499\ 530 \times 10^{-3} \pm 0.716\ 816</eq></td><td><eq>-2.362\ 255 \times 10^{-6} \pm 0.730\ 913</eq></td><td><eq>-2.666\ 849 \times 10^{-6} \pm 0.730\ 912</eq></td></tr><tr><td>500</td><td><eq>+3.029\ 965 \times 10^{-3} \pm 0.726\ 480</eq></td><td><eq>-2.668\ 530 \times 10^{-6} \pm 0.730\ 912</eq></td><td><eq>-2.668\ 473 \times 10^{-6} \pm 0.730\ 912</eq></td></tr><tr><td>1000</td><td><eq>+2.230\ 634 \times 10^{-4} \pm 0.729\ 683</eq></td><td><eq>-2.668\ 474 \times 10^{-6} \pm 0.730\ 912</eq></td><td><eq>-2.668\ 473 \times 10^{-6} \pm 0.730\ 912</eq></td></tr></table>
|
||||
1309
|
||||
09 May 2026 02:44:09
|
||||
Phys. Fluids, Vol. 16, No. 5, May 2004
|
||||
M. Sahin and R. G. Owens
|
||||
|
||||
|
||||
|
||||
FIG. 2. Reciprocal Ritz values for unbounded flow around a circular cylinder at Re546.76 computed on mesh M3 with Krylov space dimension m 5250 and shift parameter l50.5060.00i (b50.01).
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
FIG. 3. Comparison of Strouhal number versus Reynolds number for unbounded flow around a circular cylinder with other results in the literature: $( - ) .$ experimental work of Williamson ~Ref. 23!; ~¯!, numerical results of Henderson ~Ref. 25!; ~s!, numerical results of Posdziech and Grundmann ~Ref. 24!; ~h!, present ( b50.01, mesh M2!.
|
||||
|
||||
|
||||
|
||||
|
||||
FIG. 4. Comparison of drag coefficient versus Reynolds number for unbounded flow around a circular cylinder with other results in the literature: $( - ) ,$ , numerical results of Henderson ~Ref. 25!; ~s!, numerical results of Posdziech and Grundmann ~Ref. 24!; ~h!, present (b50.01, mesh M2!.
|
||||
|
||||
|
||||
# IV. FLOW PAST A CONFINED CIRCULAR CYLINDER „0.1ËbË0.9…
|
||||
Flow around a confined circular cylinder ~as opposed to the unbounded case! is an attractive benchmark problem in numerical simulation since it does not suffer from any of the difficulties associated with far-field boundary conditions ~particularly at very low Reynolds numbers! and permits the use of grid points more efficiently in smaller computational domains. Somewhat surprising, therefore, is that the only numerical linear stability analysis of Newtonian flow past a confined cylinder available in the literature would seem to be that of Chen $e t a l . ^ { 4 }$ These authors went no further than identifying the curve of neutral stability for the supercritical Hopf bifurcation at blockage ratios up to $\beta { = } 0 . 7$ . This is regrettable, because as we shall see in the paragraphs to follow, the linear stability properties of the flow become rich and therefore interesting at higher blockage ratios and Reynolds numbers than those considered by Chen et al. In the present study we consider two-dimensional flow at Reynolds numbers up to 280 and for blockage ratios in the range 0.1–0.9.
|
||||
# A. Linear stability analysis
|
||||
The curves of neutral stability computed from the GEVP with a Krylov subspace dimension $m = 2 5 0$ on mesh M2 for $\beta \in \left[ 0 . 1 , 0 . 9 \right]$ and $R e < 2 8 0$ are presented in Fig.5. Our discussion of these curves will focus on the five distinct curve sections labeled AB, BC, CD, CE, and $F G$ in the same
|
||||
1310
|
||||
09 May 2026 02:44:09
|
||||
Phys. Fluids, Vol. 16, No. 5, May 2004
|
||||
A numerical investigation of wall effects
|
||||
https://cdn-mineru.openxlab.org.cn/result/2026-05-11/e29ea199-8ede-444f-8a3a-062a2f820b92/7143dec957189df04b69e5cf9e774f69213002910deed08e6e9428081f7c6fee.jpg
|
||||
FIG. 5. Change of critical Reynolds number corresponding to both Hopf and pitchfork bifurcations with blockage ratio $\beta ,$ computed on $M 2 . A C \mathrm { : }$ Curve of neutral stability for Hopf bifurcations about symmetric solution; CD: Transition curve from asymmetric vortex shedding ~smaller $\beta )$ to a steady asymmetric solution ~larger $\beta ) ; C E \colon$ : Neutral stability curve for pitchfork bifurcation of steady symmetric solution ~smaller $\beta )$ to a steady asymmetric state ~larger $\beta ) ; F G ;$ Hopf bifurcation of an asymmetric solution ~smaller $\beta )$ to asymmetric vortex shedding ~larger $\beta ) . ~ C$ is a co-dimension 2 point where Hopf and pitchfork bifurcations occur simultaneously.
|
||||
|
||||
TABLE IV. Convergence of critical Reynolds number for different blockage ratios with $\lambda = 0 . 0 0 \pm 0 . 0 0 i$ .
|
||||
|
||||
<table><tr><td rowspan="2">Curve section(see Fig. 5)</td><td rowspan="2"><eq>\beta</eq></td><td rowspan="2">m</td><td colspan="2">M1</td><td colspan="2">M2</td><td colspan="2">M3</td><td colspan="2">Chen et al.</td></tr><tr><td><eq>Re_{\text{crit}}</eq></td><td><eq>St_{\text{crit}}</eq></td><td><eq>Re_{\text{crit}}</eq></td><td><eq>St_{\text{crit}}</eq></td><td><eq>Re_{\text{crit}}</eq></td><td><eq>St_{\text{crit}}</eq></td><td><eq>Re_{\text{crit}}</eq></td><td><eq>St_{\text{crit}}</eq></td></tr><tr><td rowspan="12">AC</td><td>0.10</td><td>500</td><td>51.00</td><td>0.1206</td><td>50.81</td><td>0.1210</td><td>50.75</td><td>0.1211</td><td>51.77</td><td>0.1116</td></tr><tr><td>0.20</td><td>250</td><td>69.86</td><td>0.1559</td><td>69.43</td><td>0.1566</td><td>69.34</td><td>0.1567</td><td>69.93</td><td>0.1559</td></tr><tr><td>0.30</td><td>250</td><td>95.24</td><td>0.2079</td><td>94.56</td><td>0.2090</td><td>94.40</td><td>0.2093</td><td>94.85</td><td>0.2085</td></tr><tr><td>0.50</td><td>250</td><td>125.23</td><td>0.3369</td><td>124.09</td><td>0.3393</td><td>123.75</td><td>0.3399</td><td>124.58</td><td>0.3382</td></tr><tr><td>0.70</td><td>250</td><td>111.32</td><td>0.4714</td><td>110.29</td><td>0.4752</td><td>110.04</td><td>0.4762</td><td>111.04</td><td>0.4744</td></tr><tr><td>0.80</td><td>250</td><td>111.45</td><td>0.5324</td><td>110.24</td><td>0.5363</td><td>109.98</td><td>0.5374</td><td></td><td></td></tr><tr><td>0.84</td><td>250</td><td>114.44</td><td>0.5530</td><td>113.69</td><td>0.5568</td><td></td><td></td><td></td><td></td></tr><tr><td>0.84</td><td>250</td><td>130.92</td><td>0.5510</td><td>126.64</td><td>0.5557</td><td></td><td></td><td></td><td></td></tr><tr><td>0.80</td><td>250</td><td>148.24</td><td>0.5324</td><td>144.19</td><td>0.5383</td><td>143.29</td><td>0.5398</td><td></td><td></td></tr><tr><td>0.76</td><td>250</td><td>169.75</td><td>0.5115</td><td>165.49</td><td>0.5186</td><td></td><td></td><td></td><td></td></tr><tr><td>0.72</td><td>250</td><td>198.94</td><td>0.4872</td><td>193.25</td><td>0.4955</td><td></td><td></td><td></td><td></td></tr><tr><td>0.70</td><td>250</td><td>218.03</td><td>0.4737</td><td>211.01</td><td>0.4827</td><td>209.40</td><td>0.4851</td><td></td><td></td></tr><tr><td rowspan="7">CE</td><td>0.70</td><td>250</td><td>221.87</td><td></td><td>216.75</td><td></td><td>215.53</td><td></td><td></td><td></td></tr><tr><td>0.72</td><td>250</td><td>210.17</td><td></td><td>205.95</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.76</td><td>250</td><td>190.65</td><td></td><td>187.01</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.80</td><td>250</td><td>173.97</td><td></td><td>169.49</td><td></td><td>168.29</td><td></td><td></td><td></td></tr><tr><td>0.84</td><td>250</td><td>161.57</td><td></td><td>158.15</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.88</td><td>250</td><td>152.93</td><td></td><td>149.84</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.90</td><td>250</td><td>147.78</td><td></td><td>145.27</td><td></td><td>144.70</td><td></td><td></td><td></td></tr><tr><td rowspan="4">CD</td><td>0.68</td><td>250</td><td>237.33</td><td>0.4596</td><td>231.06</td><td>0.4695</td><td></td><td></td><td></td><td></td></tr><tr><td>0.66</td><td>250</td><td>259.55</td><td>0.4477</td><td>253.08</td><td>0.4566</td><td></td><td></td><td></td><td></td></tr><tr><td>0.64</td><td>250</td><td>284.56</td><td>0.4351</td><td>278.01</td><td>0.4441</td><td></td><td></td><td></td><td></td></tr><tr><td>0.62</td><td>250</td><td>312.66</td><td>0.4235</td><td>306.27</td><td>0.4326</td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="6">FG</td><td>0.82</td><td>250</td><td></td><td></td><td>319.80</td><td>0.4664</td><td></td><td></td><td></td><td></td></tr><tr><td>0.82</td><td>250</td><td></td><td></td><td>227.44</td><td>0.4719</td><td></td><td></td><td></td><td></td></tr><tr><td>0.84</td><td>250</td><td>331.02</td><td>0.4954</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.84</td><td>250</td><td>214.00</td><td>0.4794</td><td>194.30</td><td>0.4979</td><td></td><td></td><td></td><td></td></tr><tr><td>0.88</td><td>250</td><td>180.43</td><td>0.5097</td><td>171.28</td><td>0.5234</td><td></td><td></td><td></td><td></td></tr><tr><td>0.90</td><td>250</td><td>169.44</td><td>0.5146</td><td>162.82</td><td>0.5202</td><td>160.50</td><td>0.5212</td><td></td><td></td></tr></table>
|
||||
1311
|
||||
09 May 2026 02:44:09
|
||||
Phys. Fluids, Vol. 16, No. 5, May 2004
|
||||
M. Sahin and R. G. Owens
|
||||
|
||||
|
||||
|
||||
FIG. 6. Change of base flow with critical Reynolds number and blockage ratio b, computed on M 2.
|
||||
|
||||
figure. The critical Reynolds numbers and corresponding Strouhal numbers ~where appropriate! for points on each curve section and computed on meshes M 1 – M 3 are supplied in Table IV.
|
||||
# 1. Curve section AB
|
||||
Validation of our numerical stability analysis and direct numerical simulations for flow past an unbounded cylinder $\beta { \approx } 0$ has been described in Sec. II. For the confined cylinder problem we have been able to compare our critical Reynolds and Strouhal numbers for the bifurcation for the symmetric state with the values for these quantities computed by Chen et $a l . ^ { 4 }$ The available results $( 0 . 1 { \leqslant } \beta { \leqslant } 0 . 7 )$ of the critical Reynolds number calculations of Chen et al. are plotted in Fig. 5 and agreement between our results and theirs over this limited section of the curve AB is excellent. Similarly excellent agreement in the computed Strouhal numbers was seen over the same range of blockage ratios, both our results and those of Chen et al. revealing a monotonic increase in the critical Strouhal number with the blockage ratio ~see Table IV!.
|
||||
Up to a blockage ratio b50.5 Table IV and the neutral stability curve AB of Fig. 5 indicate that the flow becomes more stable to two-dimensional infinitesimal disturbances as the blockage ratio increases. All along AB the flow loses stability to a Hopf bifurcation and the Strouhal number over this range of blockage ratios is increasing. Between $\beta$ 50.75 and 0.85 it may be seen from Table IV and section AB of Fig. 5 that the flow restabilizes slightly leading up to point B $( \beta = \beta _ { B } \approx 0 . 8 5 5 )$ .
|
||||
In Fig. 6 we show the streamlines of the steady base flow at seven points on the neutral stability curves. Those corresponding to point 1 are typical of those at points on and below curve AB in Fig. 5 where the steady solution is symmetric and the only recirculatory region observed is the vortex pair immediately in the wake of the cylinder itself. That is, for solutions corresponding to parameter space on and below AB in Fig. 5 no flow separation on the walls is observed.
|
||||
# 2. Curve sections BC and CE
|
||||
In Fig. 5 section BC represents the part of the critical $\beta - R e$ curve on which the time-dependent state (symmetric periodic oscillations!, passed into by crossing AB in the direction of increasing Reynolds number, restabilizes to a symmetric steady state once more. Further increases in the Reynolds number for blockage ratios in the range $\beta _ { C }$ to $\beta _ { B }$ or for choices of $\beta$ greater than $\beta _ { B }$ may result in the steady symmetric solution becoming unstable to two-dimensional perturbations via a pitchfork bifurcation into one of two asymmetric states. The curve of neutral stability for this transition is labeled C E in Fig. 5.
|
||||
|
||||
|
||||
|
||||
FIG. 7. Streamlines of unstable symmetric and stable asymmetric solutions at $R e = 1 5 0 . 0 0$ for $\beta = 0 . 9$ computed on M3.
|
||||
|
||||
1312
|
||||
09 May 2026 02:44:09
|
||||
Phys. Fluids, Vol. 16, No. 5, May 2004
|
||||
A numerical investigation of wall effects
|
||||
|
||||
|
||||
|
||||
(a)
|
||||
|
||||
|
||||
|
||||
|
||||
(b)
|
||||
|
||||
|
||||
FIG. 8. Streamlines for the disturbance velocity corresponding to ~a! the first and ~b! the second leading eigenvectors at $R e = 1 4 4 . 7 0$ for b50.9 computed on M3.
|
||||
|
||||
|
||||
The point C is a co-dimension 2 point where Hopf and pitchfork bifurcations occur simultaneously. We are able to estimate the coordinates $( \beta _ { C } , R e _ { C } )$ of this point by considering it to be the point of intersection of two straight lines drawn through the pairs of points on AC, CD, and EC that correspond to $\beta = 0 . 6 8$ and 0.7. Since b50.68 is outside the range of blockage ratios corresponding to EC the ordinate for this value of $\beta$ is computed to be that at which the leading real eigenvalue in the spectrum of linear perturbations about the ~linearly unstable! steady symmetric solution is at the origin. The critical Reynolds numbers at b50.68 and 0.7 computed on curves AC, CD, and EC are detailed in Table IV and lead to the estimate $( \beta _ { C } , R e _ { C } ) = ( 0 . 6 8 7 , 2 2 4 . 1 4 2 )$ .
|
||||
The occurrence of the transition from a symmetric steady state to an asymmetric one on CE is preceded ~in Reynolds number! by the appearance in the streamlines of a pair of downstream separation bubbles on the walls. For example, in the eigenspectrum we observed that for b50.9 and at Reynolds numbers increasing up to approximately 110 the complex conjugate pair of leading eigenvalues moved in the direction of the positive real part of the spectrum. At a Reynolds number of around 110 separation bubbles appeared on the walls and with the appearance of the separation bubbles the leading eigenpair now started to move in the opposite direction while the eigenvalue on the real axis having largest real part moved right towards the origin. In Fig. 6 we plot the streamlines on BC at points 2 and 3 to demonstrate how the size of these recirculatory regions as well as their attachment distance downstream of the cylinder increase as C is approached along the curve BC. In the context of a circular cylinder near a plane boundary such downstream separation bubbles have been observed both experimentally8 and numerically9 for cylinders sufficiently close to the boundary.
|
||||
|
||||
To gain further insight into the flow transition from steady symmetric flow ~between BC and CE and for $\beta$ $\geqslant _ { \beta _ { B } ) }$ via a pitchfork bifurcation to steady asymmetric flow ~between DE and FG) we plot in Fig. 7 the streamlines of two solutions at a Reynolds number of 150 and blockage ratio of 0.9. It may be seen from Fig. 5 that this point lies between E $( R e = 1 4 4 . 7 , \beta = \beta _ { E } = 0 . 9 )$ and $\textit { F } ( R e = 1 6 0 . 5 , \beta$ $= \beta _ { F } = 0 . 9 )$ . Thus, the symmetric solution in the upper plot in Fig. 7 is linearly unstable and the lower plot represents the streamlines of one of the stable asymmetric solutions. The disturbance velocity v in equation ~II! is, of course, solenoidal. In Fig. 8~a! we plot the streamlines associated with the disturbance field and corresponding to the dominant eigenvalue at E. The addition of a multiple of the eigenvector shown in Fig. 8~a! to the symmetric steady base flow leads to one or other of the two asymmetric steady flows, the choice dependent on the direction of circulation around the symmetric streamlines in Fig. 8~a!. An anti-clockwise direction leads to reinforcement of the lower recirculation region and reduction in the size of the upper bubble. A clockwise direction has the opposite effect. The drag coefficient associated with the steady asymmetric solution shown in Fig. 7 is slightly larger than that of the corresponding unstable symmetric one.
|
||||
1313
|
||||
09 May 2026 02:44:09
|
||||
Phys. Fluids, Vol. 16, No. 5, May 2004
|
||||
M. Sahin and R. G. Owens
|
||||
|
||||
TABLE V. Comparison with the results of Zisis and Mitsoulis ~Ref. 27! and Liu et al. ~Ref. 35! of computed total drag at Re50.0 for different blockage ratios b.
|
||||
|
||||
<table><tr><td><eq>\beta</eq></td><td>M1</td><td>M2</td><td>M3</td><td>Extrapolated</td><td>Zisis and Mitsoulis</td><td>Liu et al.</td></tr><tr><td>0.1</td><td>8.9125</td><td>8.9100</td><td>8.9089</td><td>8.9080</td><td>8.912</td><td>8.9067</td></tr><tr><td>0.2</td><td><eq>1.6221 \times 10^{1}</eq></td><td><eq>1.6215 \times 10^{1}</eq></td><td><eq>1.6211 \times 10^{1}</eq></td><td><eq>1.6200 \times 10^{1}</eq></td><td>...</td><td>...</td></tr><tr><td>0.3</td><td><eq>2.7923 \times 10^{1}</eq></td><td><eq>2.7910 \times 10^{1}</eq></td><td><eq>2.7902 \times 10^{1}</eq></td><td><eq>2.7886 \times 10^{1}</eq></td><td>...</td><td>...</td></tr><tr><td>0.5</td><td><eq>8.8354 \times 10^{1}</eq></td><td><eq>8.8294 \times 10^{1}</eq></td><td><eq>8.8263 \times 10^{1}</eq></td><td><eq>8.8227 \times 10^{1}</eq></td><td><eq>8.8207 \times 10^{1}</eq></td><td><eq>8.8227 \times 10^{1}</eq></td></tr><tr><td>0.7</td><td><eq>4.0347 \times 10^{2}</eq></td><td><eq>4.0318 \times 10^{2}</eq></td><td><eq>4.0299 \times 10^{2}</eq></td><td><eq>4.0257 \times 10^{2}</eq></td><td>...</td><td>...</td></tr><tr><td>0.9</td><td><eq>7.7057 \times 10^{3}</eq></td><td><eq>7.6988 \times 10^{3}</eq></td><td><eq>7.6959 \times 10^{3}</eq></td><td><eq>7.6941 \times 10^{3}</eq></td><td>...</td><td>...</td></tr></table>
|
||||
# 3. Curve section CD
|
||||
The section of the neutral stability curves labeled CD in Fig. 5 represents a transition curve with increasing b from periodic vortex shedding to the left of this curve ~smaller b! to a steady asymmetric state ~larger b!. At point C the steady solution is symmetric but moving along the curve CD towards D causes the growth of one of the recirculatory wall regions relative to the other.
|
||||
Since all curves of neutral stability in Fig. 5 have been determined using a linear stability analysis about a steady flow as described in Sec. II we have been unable to plot the precise boundaries of the transition region that must exist from symmetric oscillations to asymmetric oscillations as the curve CD is approached in parameter space from the left ~smaller b!.
|
||||
# 4. Curve section FG
|
||||
Finally, the steady asymmetric solution of the region between curves DE and FG can become unstable via a Hopf bifurcation to asymmetric vortex shedding ~see discussion in Sec. IV B!. The transition curve is plotted as FG and the streamlines of two steady base flows at points 6 and 7 of this curve are shown in Fig. 6. In Fig. 8~b! we plot the streamlines of the disturbance velocity corresponding to the leading complex eigenvalue pair in the spectrum at E. It is the addition of a mode similar in form to this ~but at a higher Reynolds number! that leads to vortex shedding about the asymmetric state. If the Reynolds number is further increased on the curve FG additional separation bubbles appear on the wall further downstream. We also remark that on the curve FG the separation bubble just behind the cylinder is generally shorter and more rounded than that computed on the curves BC and CD.
|
||||
A strong parallel is thus seen in the present results with those of numerous other authors ~see, for example, those of Battaglia et al., 28 Drikakis,29 Fearn et al., 30 Hawa and Rusak,31 Mishra and Jayaraman,32 and Oliveira33! for flows through both two-dimensional and three-dimensional symmetric expansions. All the cited authors report that steady flow with symmetric recirculatory regions through an expansion geometry encounters a supercritical pitchfork bifurcation at a certain Reynolds number ~dependent, of course, on the channel geometry! and becomes asymmetric. The difference in the streamwise attachment length of the two recirculatory regions ~still in the steady regime! becomes larger as the Reynolds number is further increased from the critical value. In the two-dimensional case, increasing the expansion ratio decreases the critical Reynolds number.28 In the threedimensional case, Schrek and Scha¨fer34 found that fixing the expansion ratio at 1:3 and decreasing the width of the channel relative to the downstream channel height from ` ~twodimensional flow! through 5 to 2 resulted in a stabilization of the flow.
|
||||
|
||||
It will be noted from Table IV that for the choices of β=0.82 and 0.84 critical Reynolds numbers of 319.8 and 331.02, respectively, are added to those that are shown in Fig. 5. This is to indicate how the curve FG would continue if the range of Reynolds numbers were to be extended in Fig. 5, although at these higher Reynolds numbers it is highly unlikely that the flow would in reality remain twodimensional.
|
||||
# B. Direct numerical simulations
|
||||
A few verifications were performed on the results of our direct numerical simulations in order to establish their reliability. First, Strouhal numbers near the critical Reynolds numbers corresponding to the onset of periodic vortex shedding and computed with direct numerical simulation were found to be in good agreement in a couple of cases with those predicted on the basis of the eigenvalue analysis of
|
||||
|
||||
|
||||
|
||||
FIG. 9. Computed drag coefficient versus Reynolds number at blockage ratios b5 0.1, 0.3, 0.5, 0.7, and 0.9.
|
||||
|
||||
1314
|
||||
09 May 2026 02:44:09
|
||||
Phys. Fluids, Vol. 16, No. 5, May 2004
|
||||
A numerical investigation of wall effects
|
||||
|
||||
|
||||
|
||||
FIG. 10. Change of time-dependent flow with Reynolds number and blockage ratio b, computed on M 2.
|
||||
|
||||
Sec. IV A. For example, at a blockage ratio of $\beta = 0 . 3$ and at $R e = 1 0 0$ the corresponding Strouhal number was computed from the lift coefficient data over extended time intervals and found to be equal to 0.2115. This compares well with the value of 0.2090 supplied in Table IV and computed at $R e _ { \mathrm { c r i t } } { = } 9 4 . 5 6$ . Second, we present in Table V results of computations using the three meshes M 1 to M 3 of the drag on the cylinder for various blockage ratios. These are compared with the recent numerical data of Zisis and Mitsoulis27 and Liu $e t a l . ^ { 3 5 }$ and the agreement is convincing. It should be further added that the drag result of Liu et al. for $\beta { = } 0 . 1$ is within 0.007% of the theoretically predicted value of Faxe´n.36 The drag coefficient versus Reynolds number is given in Fig. 9 for several blockage ratios. At all the blockage ratios considered here the drag coefficient behaves like $1 / R e$ at low Reynolds numbers. As the blockage ratio increases the range of values of Re over which this remains true gets smaller.
|
||||
|
||||
|
||||
|
||||
|
||||
(a)t=0
|
||||
|
||||
|
||||
|
||||
|
||||
(b)t=T/3
|
||||
|
||||
|
||||
|
||||
|
||||
(c)t=2T/3
|
||||
|
||||
|
||||
FIG. 11. Vorticity contours of the periodic flow at $R e = 2 0 0 . 0 0$ and b50.5 ~point 1 in Fig. 10! computed on M2. t50 corresponds to the solution having minimum lift coefficient and the period T'2.85.
|
||||
|
||||
1315
|
||||
09 May 2026 02:44:09
|
||||
Phys. Fluids, Vol. 16, No. 5, May 2004
|
||||
M. Sahin and R. G. Owens
|
||||
|
||||
|
||||
|
||||
(a)
|
||||
|
||||
|
||||
|
||||
|
||||
(b)
|
||||
|
||||
|
||||
|
||||
|
||||
（c）
|
||||
|
||||
|
||||
|
||||
|
||||
(d）
|
||||
|
||||
|
||||
|
||||
|
||||
(e）
|
||||
|
||||
|
||||
FIG. 12. Phase space plots of lift and drag coefficients parametrized with nondimensional time $( t U _ { \mathrm { m a x } } / D )$ , computed on M2. ~a! $R e = 2 0 0 . 0 0$ and $\beta { = } 0 . 5$ ~point 1 in Fig. 10!, ~b! $R e = 2 0 0 . 0 0$ and $\beta { = } 0 . 7$ ~point 2 in Fig. 10!, ~c! $R e = 2 0 0 . 0 0$ and $\beta { = } 0 . 9$ ~point 3 in Fig. 10!, ~d! $R e = 2 0 0 . 0 0$ and $\beta { = } 0 . 8$ ~point 4 in Fig. 10!, ~e! $R e = 1 6 0 . 0 0$ and $\beta { = } 0 . 8$ ~point 5 in Fig. 10!.
|
||||
|
||||
1316
|
||||
09 May 2026 02:44:09
|
||||
Phys. Fluids, Vol. 16, No. 5, May 2004
|
||||
A numerical investigation of wall effects
|
||||
|
||||
|
||||
|
||||
(a)t=0
|
||||
|
||||
|
||||
|
||||
|
||||
(b) $t = T / 3$
|
||||
|
||||
|
||||
|
||||
|
||||
$( \mathrm { c } ) \ t = 2 T / 3$
|
||||
|
||||
|
||||
FIG. 13. Vorticity contours of the periodic flow at $R e = 2 0 0 . 0 0$ and $\beta { = } 0 . 7$ ~point 2 in Fig. 10! computed on $M 2 . \ t = 0$ corresponds to the solution having minimum lift coefficient and the period T'2.05.
|
||||
|
||||
|
||||
|
||||
|
||||
(a)t=0
|
||||
|
||||
|
||||
|
||||
|
||||
$( \mathbf { b } ) \ t = T / 3$
|
||||
|
||||
|
||||
|
||||
|
||||
$( \mathrm { c } ) \ t = 2 T / 3$
|
||||
|
||||
|
||||
FIG. 14. Vorticity contours of the periodic solution at $R e = 2 0 0 . 0 0$ and $\beta { = } 0 . 9$ ~point 3 in Fig. 10! computed on $M 2 , t = 0$ corresponds to the solution having minimum lift coefficient and the period T'1.88.
|
||||
|
||||
|
||||
|
||||
|
||||
(a） $t = 0$
|
||||
|
||||
|
||||
|
||||
|
||||
(b) $t = T / 3$
|
||||
|
||||
|
||||
|
||||
|
||||
$( \mathrm { c } ) \ t = 2 T / 3$
|
||||
|
||||
|
||||
FIG. 15. Vorticity contours of the periodic solution at $R e = 2 0 0 . 0 0$ and $\beta { = } 0 . 8$ ~point 4 in Fig. 10! computed on $M 2 , t = 0$ corresponds to the solution having minimum lift coefficient and the period $T { \approx } 1 . 8 1 5 .$ .
|
||||
|
||||
1317
|
||||
09 May 2026 02:44:09
|
||||
Phys. Fluids, Vol. 16, No. 5, May 2004
|
||||
M. Sahin and R. G. Owens
|
||||
|
||||
|
||||
|
||||
(a）t=0
|
||||
|
||||
|
||||
|
||||
|
||||
(b)t=T/3
|
||||
|
||||
|
||||
|
||||
|
||||
(c)t=2T/3
|
||||
|
||||
|
||||
FIG. 16. Vorticity contours of the periodic solution at $R e = 1 6 0 . 0 0$ and b50.8 ~point 5 in Fig. 10! computed on M2. t50 corresponds to the solution having minimum lift coefficient and the period T'1.806.
|
||||
|
||||
In order to elucidate the variation in the critical Reynolds number with blockage ratio observed from the eigenvalue analysis of Sec. IV A we used direct numerical simulation to investigate the wake structure at five different locations in the $\beta - R e$ parameter space and labeled 1–5 in Fig. 10. As an aside, and before discussing our results in detail, we simply note that at blockage ratios $\beta { < } 0 . 5$ direct numerical simulations revealed that the vortex shedding over the cylinder was quite similar to that of the unbounded case, although the vortex street is shorter due to shear in the free stream. We also note that although blockage effects are expected to delay transition of the cylinder wake to threedimensional flow, it is possible that at some of the points labeled 1–4 in Fig. 10 the local velocity is so high that a three-dimensional transition occurs for the highly decelerated, separated boundary layers on the channel walls. Verification of this will have to await fully three-dimensional simulations, however.
|
||||
Time-dependent solutions are presented in Fig. 11 for $\beta = 0 . 5$ at a Reynolds number of 200 (point 1 of Fig.10).At this Reynolds number the flow has lost its stability to twodimensional disturbances and has become time-periodic with ${ S t = 0 . 3 5 1 3 }$ , which is higher than for unbounded flow around a circular cylinder $( S t = 0 . 1 9 7 7 )$ . The sequences of three ‘‘snapshots’’ in Fig. 11 are taken at times t50, T/3 and 2T/3, the nondimensional period T being approximately equal to 2.85 and determined from the lift coefficient data over long time periods. t50 corresponds to a minimum in the lift coefficient once fully periodic vortex shedding is established. It may be seen from Figs. 11~a!–11~c! that vortex shedding occurs both from the cylinder and the channel walls. As these vortices move downstream the trajectories of clockwise vortices shed from the upper part of the cylinder cross those having opposite sign ~and shed from the lower part of the cylinder! so that wall proximity effects are seen to give rise to a reverse von Ka´rma´n street. The same phenomenon has been documented by other authors.9 In Fig. 12~a! we show the $C _ { d } { - } C _ { l }$ phase space plot at point 1 of Fig. 10, once fully periodic conditions have been established. The average lift is zero and $C _ { d }$ and $C _ { l }$ are both symmetric in the rising and falling parts of each cycle.
|
||||
|
||||
At a blockage ratio of 0.7 the flow is periodic at a Reynolds number of 200.00 ~point 2 of Fig. 10! with St 50.4881. The vorticity contours at times t50, T/3 and 2T/3 are shown in Figs. 13~a!–13~c!, with the period T'2.05. Unlike in the case of $\beta = 0 . 5$ described in the paragraph above, vortex shedding from the cylinder seems to be almost suppressed at this Reynolds number, due to the proximity of point 2 to the curve of neutral stability BC. However, there are very weak vortices shed from both upper and lower lateral walls. These are well separated from each other and their interaction is weak. The phase space plot of the lift and drag coefficients at this blockage ratio are shown in Fig. 12~b!. Although the time-averaged value of the drag coefficient $C _ { d }$ has increased it is notable that the amplitude of the $C _ { d }$ oscillations is an order of magnitude less than that seen at $\beta$ 50.5.
|
||||
At a blockage ratio of 0.9 and $R e > 1 6 0 . 5$ the flow is unsteady and very strong vortices are shed from both the cylinder and the walls. The computed streamlines and vorticity contours at $R e = 2 0 0 . 0 0$ ~point 3 of Fig. 10! are shown in Fig. 14, each separated from the previous in the series by a third of a period T/3. At this Reynolds number the flow is periodic with St50.5314 (T'1.88) and this compares reasonably well ~see Table IV! with the Strouhal number of 0.5202 computed from the GEVP on mesh M 2 at the critical Reynolds number $R e _ { \mathrm { c r i t } } { = } 1 6 2 . 8 2$ . Vortices having the same sign merge just behind the cylinder and are then transported downstream. However, the vortex street formed behind the cylinder is quite different from the well-known von Ka´rma´n street in that very strong opposite-sign vortices with smaller structure move downstream and interact with the wall, creating strong vortices there. Although the streamlines of the eigenvector shown in Fig. 8~b! will be slightly modified by the time the Reynolds number reaches 200, the size of the cellular structures in the wake of the eigenvector are close to what is seen in the direct numerical solution for the vorticity in Fig. 14. Oscillations are now about an asymmetric state and lead to drastic increases in both the lift and drag coefficient values. Their variation over a cycle with dimensionless time is shown in Fig. 12~c!. In this figure it may be seen that not only is the time-averaged value of $C _ { d }$ greater than for the two previous blockage ratios considered, but the amplitudes of oscillation of both coefficients has dramatically increased after a tendency observed up to $\beta { \approx } 0 . 7 5$ of successively diminishing amplitudes. We also note for the first time that the lift coefficient is no longer symmetric in the rising and falling parts of each cycle ~the ‘‘figure of $8 ^ { \circ }$ is distorted!, although the average value of $C _ { l }$ over one cycle is zero.
|
||||
1318
|
||||
09 May 2026 02:44:09
|
||||
Phys. Fluids, Vol. 16, No. 5, May 2004
|
||||
A numerical investigation of wall effects
|
||||
|
||||
The flow behavior in the parametric region between curves BD and FG is particularly interesting since although the steady base flow solutions are linearly stable they appear to be finitely unstable. Although our primary concern in this paper is the study of wall effects on the linear stability of flow past a confined cylinder, in Fig. 10 we plot the instantaneous streamlines of two pairs of different solutions, both pairs being for a blockage ratio $\beta { = } 0 . 8$ but corresponding to two different values of $R e \colon$ at point 4 $R e = 2 0 0$ whereas at point 5 $R e = 1 6 0$ . Shown in the uppermost plots at points 4 and 5 are the steady and linearly stable solutions, but these may be sent into permanently unsteady states by running, for example, the time-dependent code with the geometrically rescaled steady base flow of point $2 \ ( \beta = 0 . 7 , R e = 2 0 0 )$ as the initial ‘‘guess.’’ Snapshots of unsteady flow at times t $= 0 , T / 3$ , and 2T/3 at point 4 are presented in Figs. 15~a!– 15~c!, with the period $T { \approx } 1 . 8 1 5$ . Unlike the other unsteady cases considered so far in this paper, the recirculation regions on the upper wall are much larger than those on the lower wall and they can move far downstream. Additionally, the recirculatory region at the lower part of the cylinder is larger than that at the upper part. This asymmetry may be seen from the phase space plot in Fig. 12~d!. The time average of the lift coefficient is no longer zero and the lift curve is asymmetric in the rising and falling part of each cycle. However, if the Reynolds number is chosen equal to 160 ~point 5 in Fig. 10! the flow becomes symmetric again since point 5 is below the curve CE. The computed vorticity contours are given in Fig. 16 with T'1.806. At this point the flow structure is quite similar to that of point 2 with vortex shedding from the upper and lower walls. The lift and drag coefficients which are supplied in the phase space plot in Fig. 12~e! are also similar to those of point 2 with zero time average of the lift coefficient.
|
||||
Numerical experiments at a blockage ratio of 0.9 and a Reynolds number of 500 indicated that the flow had become chaotic. However, it seems unlikely that the flow is still twodimensional at this Reynolds number and presentation of our results will have to await a fully three-dimensional analysis.
|
||||
# V. CONCLUSIONS
|
||||
In this paper we have computed with greater accuracy and over a larger range of blockage ratios than has proved possible in the past the effects on the drag and linear stability of lateral wall proximity for flow past a cylinder at Reynolds numbers up to 280.
|
||||
Some of the rich and complex dynamics of the system for sufficiently high Reynolds numbers and blockage ratios have been uncovered and discussed. In particular, we have found that for $R e \leqslant 2 8 0$ and $\beta { \leqslant } 0 . 9$ there are ~at least! three separate curves of neutral stability: ~a! Hopf bifurcation of a symmetric state, ~b! pitchfork bifurcation of a symmetric state to one of two asymmetric states, and ~c! Hopf bifurcation of an asymmetric state leading to asymmetric oscillations thereafter. In addition, we have drawn attention to a transition region from symmetric vortex shedding to asymmetric vortex shedding with increasing blockage ratio. Further increases in the blockage ratio ~crossing CD in Fig. 10! leads to restabilization to a steady asymmetric solution.
|
||||
A co-dimension 2 point where pitchfork and Hopf bifurcations occur simultaneously has been identified and a region in parameter space ~either side of the locus of the pitchfork bifurcation! seems to exist where the steady solution is linearly stable but unstable to finite two-dimensional perturbations.
|
||||
# ACKNOWLEDGMENTS
|
||||
The authors wish to thank Peter Monkewitz for sharing with them interesting and illuminating insights into the primary and secondary instability mechanisms. The work of the first author was supported by the Swiss National Science Foundation, Grant No. 21-61865.00.
|
||||
|
||||
|
||||
1C. H. K. Williamson, ‘‘Vortex dynamics in the cylinder wake,’’ Annu. Rev. Fluid Mech. 28, 477 ~1996!.
|
||||
|
||||
|
||||
|
||||
|
||||
2P. Anagnostopoulos, G. Iliadis, and S. Richardson, ‘‘Numerical study of the blockage effect on viscous flow past a circular cylinder,’’ Int. J. Numer. Methods Fluids 22, 1061 ~1996!.
|
||||
|
||||
|
||||
|
||||
|
||||
3M. Coutanceau and R. Bouard, ‘‘Experimental determination of the main features of the viscous flow in the wake of a circular cylinder in uniform translation. Part 1. Steady flow,’’ J. Fluid Mech. 79, 231 ~1977!.
|
||||
|
||||
|
||||
|
||||
|
||||
4J.-H. Chen, W. G. Pritchard, and S. J. Tavener, ‘‘Bifurcation for flow past a cylinder between parallel planes,’’ J. Fluid Mech. 284, 23 ~1995!.
|
||||
|
||||
|
||||
|
||||
|
||||
5M. Behr, S. Hastreiter, S. Mittal, and T. E. Tezduyar, ‘‘Incompressible flow past a circular cylinder: Dependence of the computed flow field on the location of the lateral boundaries,’’ Comput. Methods Appl. Mech. Eng. 123, 309 ~1995!.
|
||||
|
||||
|
||||
|
||||
|
||||
6P. K. Stansby and A. Slaouti, ‘‘Simulation of vortex shedding including blockage by the random-vortex and other methods,’’ Int. J. Numer. Methods Fluids 17, 1003 ~1993!.
|
||||
|
||||
|
||||
|
||||
|
||||
7C. Lei, L. Cheng, and K. Kavanagh, ‘‘Re-examination of the effect of a plane boundary on force and vortex shedding of a circular cylinder,’’ J. Wind. Eng. Ind. Aerodyn. 80, 263 ~1999!.
|
||||
|
||||
|
||||
|
||||
|
||||
8P. W. Bearman and M. M. Zdravkovich, ‘‘Flow around a circular cylinder near a plane boundary,’’ J. Fluid Mech. 89, 33 ~1978!.
|
||||
|
||||
|
||||
|
||||
|
||||
9L. Zovatto and G. Pedrizzetti, ‘‘Flow around a circular cylinder between parallel walls,’’ J. Fluid Mech. 440, 1 ~2001!.
|
||||
|
||||
|
||||
|
||||
|
||||
10C. P. Jackson, ‘‘A finite-element study of the onset of vortex shedding in flow past variously shaped bodies,’’ J. Fluid Mech. 182, 23 ~1987!.
|
||||
|
||||
|
||||
|
||||
|
||||
11M. Sahin, ‘‘Solution of the incompressible unsteady Navier–Stokes equations only in terms of the velocity components,’’ Int. J. Comput. Fluid Dyn. 17, 199 ~2003!.
|
||||
|
||||
|
||||
|
||||
|
||||
12M. Sahin and R. G. Owens, ‘‘A novel fully-implicit finite volume method applied to the lid-driven cavity problem. Part I. High Reynolds number
|
||||
|
||||
|
||||
1319
|
||||
09 May 2026 02:44:09
|
||||
Phys. Fluids, Vol. 16, No. 5, May 2004
|
||||
M. Sahin and R. G. Owens
|
||||
|
||||
|
||||
flow calculations,’’ Int. J. Numer. Methods Fluids 42, 57 ~2003!.
|
||||
|
||||
|
||||
|
||||
|
||||
13M. Sahin and R. G. Owens, ‘‘A novel fully-implicit finite volume method applied to the lid-driven cavity problem. Part II. Linear stability analysis,’’ Int. J. Numer. Methods Fluids 42, 79 ~2003!.
|
||||
|
||||
|
||||
|
||||
|
||||
14A. Kourta, M. Braza, P. Chassaing, and H. Haminh, ‘‘Numerical analysis of a natural and excited two-dimensional mixing layer,’’ AIAA J. 25, 279 ~1987!.
|
||||
|
||||
|
||||
|
||||
|
||||
15G. Jin and M. Braza, ‘‘A nonreflecting outlet boundary condition for incompressible unsteady Navier–Stokes calculations,’’ J. Comput. Phys. 107, 239 ~1993!.
|
||||
|
||||
|
||||
|
||||
|
||||
16W. E. Arnoldi, ‘‘The principle of minimized iterations in the solution of the matrix eigenvalue problem,’’ Q. Appl. Math. 9, 17 ~1951!.
|
||||
|
||||
|
||||
|
||||
|
||||
17Y. Saad, ‘‘Variations on Arnoldi’s method for computing eigen elements of large unsymmetric matrices,’’ Linear Algebr. Appl. 34, 269 ~1980!.
|
||||
|
||||
|
||||
|
||||
|
||||
18R. Natarajan, ‘‘An Arnoldi-based iterative scheme for nonsymmetric matrix pencils arising in finite element stability problems,’’ J. Comput. Phys. 100, 128 ~1992!.
|
||||
|
||||
|
||||
|
||||
|
||||
19P. R. Amestoy, I. S. Duff, and J.-Y. L’Excellent, ‘‘Multifrontal parallel distributed symmetric and unsymmetric solvers,’’ Comput. Methods Appl. Mech. Eng. 184, 501 ~2000!.
|
||||
|
||||
|
||||
|
||||
|
||||
20P. R. Amestoy, I. S. Duff, J. Koster, and J.-Y. L’Excellent, ‘‘A fully asynchronous multifrontal solver using distributed dynamic scheduling,’’ SIAM J. Matrix Anal. Appl. 23, 15 ~2001!.
|
||||
|
||||
|
||||
|
||||
|
||||
21J. L. Steger and R. L. Sorenson, ‘‘Automatic mesh-point clustering near a boundary in grid generation with elliptic partial differential equations,’’ J. Comput. Phys. 33, 405 ~1979!.
|
||||
|
||||
|
||||
|
||||
|
||||
22Y. Ding and M. Kawahara, ‘‘Three-dimensional linear stability analysis of incompressible viscous flows using the finite element method,’’ Int. J. Numer. Methods Fluids 31, 451 ~1999!.
|
||||
|
||||
|
||||
|
||||
|
||||
23C. H. K. Williamson, ‘‘Oblique and parallel modes of vortex shedding in the wake of a circular cylinder at low Reynolds numbers,’’ J. Fluid Mech. 206, 579 ~1989!.
|
||||
|
||||
|
||||
|
||||
|
||||
24O. Posdziech and R. Grundmann, ‘‘Numerical simulation of the flow around an infinitely long circular cylinder in the transition regime,’’ Theor. Comput. Fluid Dyn. 15, 121 ~2001!.
|
||||
|
||||
|
||||
|
||||
|
||||
25R. D. Henderson, ‘‘Detail of the drag curve near the onset of vortex shedding,’’ Phys. Fluids 7, 2102 ~1995!.
|
||||
|
||||
|
||||
|
||||
|
||||
26B. Fornberg, ‘‘A numerical study of steady viscous flow past a circular cylinder,’’ J. Fluid Mech. 98, 819 ~1980!.
|
||||
|
||||
|
||||
|
||||
|
||||
27Th. Zisis and E. Mitsoulis, ‘‘Viscoplastic flow around a cylinder kept between parallel plates,’’ J. Non-Newtonian Fluid Mech. 105, 1 ~2002!.
|
||||
|
||||
|
||||
|
||||
|
||||
28F. Battaglia, S. J. Tavener, A. K. Kulkarni, and C. L. Merkle, ‘‘Bifurcation of low Reynolds number flows in symmetric channels,’’ AIAA J. 35, 99 ~1997!.
|
||||
|
||||
|
||||
|
||||
|
||||
29D. Drikakis, ‘‘Bifurcation phenomena in incompressible sudden expansion flows,’’ Phys. Fluids 9, 76 ~1997!.
|
||||
|
||||
|
||||
|
||||
|
||||
30R. M. Fearn, T. Mullin, and K. A. Cliffe, ‘‘Nonlinear flow phenomena in a symmetric sudden expansion,’’ J. Fluid Mech. 211, 595 ~1990!.
|
||||
|
||||
|
||||
|
||||
|
||||
31T. Hawa and Z. Rusak, ‘‘The dynamics of a laminar flow in a symmetric channel with a sudden expansion,’’ J. Fluid Mech. 436, 283 ~2001!.
|
||||
|
||||
|
||||
|
||||
|
||||
32S. Mishra and K. Jayaraman, ‘‘Asymmetric flows in planar symmetric channels with large expansion ratio,’’ Int. J. Numer. Methods Fluids 38, 945 ~2002!.
|
||||
|
||||
|
||||
|
||||
|
||||
33P. J. Oliveira, ‘‘Asymmetric flows of viscoelastic fluids in symmetric planar expansion geometries,’’ J. Non-Newtonian Fluid Mech. 114, 33 ~2003!.
|
||||
|
||||
|
||||
|
||||
|
||||
34E. Schreck and M. Scha¨fer, ‘‘Numerical study of bifurcation in threedimensional sudden channel expansions,’’ Comput. Fluids 29, 583 ~2000!.
|
||||
|
||||
|
||||
|
||||
|
||||
35A. W. Liu, D. E. Bornside, R. C. Armstrong, and R. A. Brown, ‘‘Viscoelastic flow of polymer solutions around a periodic, linear array of cylinders: comparisons of predictions for microstructure and flow fields,’’ J. Non-Newtonian Fluid Mech. 77, 153 ~1998!.
|
||||
|
||||
|
||||
|
||||
|
||||
36O. H. Faxe´n, ‘‘Forces exerted on a rigid cylinder in a viscous fluid between two parallel fixed planes,’’ R. Swed. Acad. Eng. Sci. 187, 1 ~1946!.
|
||||
|
||||
|
||||
1320
|
||||
File diff suppressed because it is too large
Load Diff
@@ -1,755 +0,0 @@
|
||||
#!/usr/bin/env python3
|
||||
"""
|
||||
Stability Matrix Test
|
||||
=====================
|
||||
Tests three collision models (SRT/TRT/MRT) at low and high Re (with/without LES),
|
||||
plus Esoteric-Pull streaming at low Re with SRT.
|
||||
|
||||
Outputs:
|
||||
- Flow-field images (velocity, vorticity, streamlines) for each case
|
||||
- Diagnostic JSON with stability metrics
|
||||
- EsoPull vs double-buffer comparison plots
|
||||
|
||||
Usage:
|
||||
python3 tests/test_stability_matrix.py [--device 0] [--steps 2000]
|
||||
"""
|
||||
|
||||
import argparse
|
||||
import json
|
||||
import math
|
||||
import os
|
||||
import struct
|
||||
import sys
|
||||
import time
|
||||
|
||||
sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "src"))
|
||||
|
||||
import matplotlib
|
||||
matplotlib.use("Agg")
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
import pycuda.driver as cuda
|
||||
|
||||
from CelerisLab.cuda import compiler
|
||||
|
||||
# ---------------------------------------------------------------------------
|
||||
# Constants
|
||||
# ---------------------------------------------------------------------------
|
||||
FLUID = 0x01
|
||||
SOLID = 0x02
|
||||
OBSTACLE = 0x20 # fixed: was 0x04
|
||||
|
||||
COLLISION_NAMES = {0: "SRT", 1: "TRT", 2: "MRT"}
|
||||
|
||||
|
||||
# ---------------------------------------------------------------------------
|
||||
# Helpers
|
||||
# ---------------------------------------------------------------------------
|
||||
def compute_vis_omega(re, diameter, u0):
|
||||
vis = u0 * diameter / re
|
||||
omega = 1.0 / (3.0 * vis + 0.5)
|
||||
return vis, omega
|
||||
|
||||
|
||||
def lattice_weights(nq):
|
||||
if nq == 9:
|
||||
return np.array([4/9] + [1/9]*4 + [1/36]*4, dtype=np.float32)
|
||||
if nq == 19:
|
||||
return np.array([1/3] + [1/18]*6 + [1/36]*12, dtype=np.float32)
|
||||
raise ValueError(f"nq={nq}")
|
||||
|
||||
|
||||
def build_flags_2d(nx, ny, cx, cy, radius):
|
||||
flag = np.ones(nx * ny, dtype=np.uint8) * FLUID
|
||||
for y in range(ny):
|
||||
for x in range(nx):
|
||||
k = y * nx + x
|
||||
if y == 0 or y == ny - 1 or x == 0 or x == nx - 1:
|
||||
flag[k] = SOLID
|
||||
elif (x - cx)**2 + (y - cy)**2 < radius**2:
|
||||
flag[k] = OBSTACLE
|
||||
return flag
|
||||
|
||||
|
||||
def set_macros(nx, ny, dim, nq, vis, u0, collision_model, use_les, streaming_model,
|
||||
omega_collision_max=1.999, inlet_profile=1, trt_magic_param=0.1875,
|
||||
les_cs=0.16):
|
||||
"""Write config/*.h files used by kernel_v2.cu."""
|
||||
cfg_dir = os.path.join(os.path.dirname(compiler.kernel_path("config.h")), "config")
|
||||
|
||||
# config_grid.h
|
||||
with open(os.path.join(cfg_dir, "config_grid.h"), "w") as f:
|
||||
f.write(f"""\
|
||||
// AUTO-GENERATED by test_stability_matrix.py
|
||||
#ifndef CELERIS_CONFIG_GRID_H
|
||||
#define CELERIS_CONFIG_GRID_H
|
||||
#define NT 128
|
||||
#define MULT_GPU 0
|
||||
#define NX {nx}
|
||||
#define NY {ny}
|
||||
#define NZ 1
|
||||
#define DIM {dim}
|
||||
#define NQ {nq}
|
||||
#endif
|
||||
""")
|
||||
|
||||
# config_physics.h
|
||||
with open(os.path.join(cfg_dir, "config_physics.h"), "w") as f:
|
||||
f.write(f"""\
|
||||
// AUTO-GENERATED by test_stability_matrix.py
|
||||
#ifndef CELERIS_CONFIG_PHYSICS_H
|
||||
#define CELERIS_CONFIG_PHYSICS_H
|
||||
#define LBtype float
|
||||
#define VIS {vis:.10f}
|
||||
#define RHO 1.0
|
||||
#define U0 {u0}
|
||||
#define PI 3.141592653589793238
|
||||
#define FLUID 0x01
|
||||
#define SOLID 0x02
|
||||
#define GAS 0x04
|
||||
#define INTERFACE 0x08
|
||||
#define SENSOR 0x10
|
||||
#define OBSTACLE 0x20
|
||||
#define V_TAYLOR 1
|
||||
#endif
|
||||
""")
|
||||
|
||||
# config_method.h
|
||||
with open(os.path.join(cfg_dir, "config_method.h"), "w") as f:
|
||||
f.write(f"""\
|
||||
// AUTO-GENERATED by test_stability_matrix.py
|
||||
#ifndef CELERIS_CONFIG_METHOD_H
|
||||
#define CELERIS_CONFIG_METHOD_H
|
||||
#define COLLISION_MODEL {collision_model}
|
||||
#define STREAMING_MODEL {streaming_model}
|
||||
#define STORE_PRECISION 0
|
||||
#define USE_DDF_SHIFTING 0
|
||||
#define USE_LES {int(use_les)}
|
||||
#define LES_CS {les_cs:.6f}f
|
||||
#define INLET_PROFILE {int(inlet_profile)}
|
||||
#define OUTLET_MODE 0
|
||||
#define OUTLET_BLEND_ALPHA 0.700f
|
||||
#define OUTLET_BACKFLOW_CLAMP 1
|
||||
#define OMEGA_COLLISION_MIN 0.01f
|
||||
#define OMEGA_COLLISION_MAX {float(omega_collision_max):.3f}f
|
||||
#define TRT_MAGIC_PARAM {float(trt_magic_param):.6f}f
|
||||
#endif
|
||||
""")
|
||||
|
||||
# config_objects.h
|
||||
with open(os.path.join(cfg_dir, "config_objects.h"), "w") as f:
|
||||
f.write("""\
|
||||
// AUTO-GENERATED by test_stability_matrix.py
|
||||
#ifndef CELERIS_CONFIG_OBJECTS_H
|
||||
#define CELERIS_CONFIG_OBJECTS_H
|
||||
#define N_OBJS 0
|
||||
#endif
|
||||
""")
|
||||
|
||||
|
||||
def pack_d_params(nx, ny, omega, u0):
|
||||
"""Pack LBMParams struct for __constant__ memory upload."""
|
||||
return struct.pack(
|
||||
"IIIQfffffffI",
|
||||
nx, ny, 1, # Nx, Ny, Nz
|
||||
nx * ny, # N
|
||||
omega, # omega
|
||||
1.1, # omega_bulk
|
||||
0.0, 0.0, 0.0, # fx, fy, fz
|
||||
1.0, # rho_ref
|
||||
u0, # u_inlet
|
||||
0, # n_objects
|
||||
)
|
||||
|
||||
|
||||
def impose_rest_on_nonfluid(flag, host_ddf, nq, nx, ny):
|
||||
w = lattice_weights(nq)
|
||||
f = host_ddf.reshape(nq, ny, nx)
|
||||
nonfluid = flag.reshape(ny, nx) != FLUID
|
||||
for i in range(nq):
|
||||
f[i, nonfluid] = w[i]
|
||||
return host_ddf
|
||||
|
||||
|
||||
def compute_macros_2d(host_ddf, nq, nx, ny, flag):
|
||||
"""Compute rho, ux, uy from DDF."""
|
||||
cx9 = [0, 1, -1, 0, 0, 1, -1, 1, -1]
|
||||
cy9 = [0, 0, 0, 1, -1, 1, -1, -1, 1]
|
||||
f = host_ddf.reshape(nq, ny, nx)
|
||||
rho = np.sum(f, axis=0)
|
||||
ux = np.zeros_like(rho)
|
||||
uy = np.zeros_like(rho)
|
||||
for i in range(nq):
|
||||
ux += cx9[i] * f[i]
|
||||
uy += cy9[i] * f[i]
|
||||
rho_safe = np.where(np.abs(rho) > 1e-12, rho, 1.0)
|
||||
ux /= rho_safe
|
||||
uy /= rho_safe
|
||||
return rho, ux, uy
|
||||
|
||||
|
||||
def diagnose(rho, ux, uy, flag, nx, ny):
|
||||
"""Compute stability diagnostics."""
|
||||
fluid = flag.reshape(ny, nx) == FLUID
|
||||
nan_count = int(np.isnan(rho).sum())
|
||||
rho_min = float(np.nanmin(rho))
|
||||
rho_max = float(np.nanmax(rho))
|
||||
mass = float(np.nansum(rho[fluid]))
|
||||
vel = np.sqrt(ux**2 + uy**2)
|
||||
|
||||
# Ma check
|
||||
ma_max = float(np.nanmax(vel[fluid])) * math.sqrt(3.0) if np.any(fluid) else 0.0
|
||||
|
||||
# Vorticity RMS in wake region
|
||||
vort = np.gradient(uy, axis=1) - np.gradient(ux, axis=0)
|
||||
wake_mask = fluid & (np.arange(nx)[None, :] > nx // 3)
|
||||
vort_rms = float(np.sqrt(np.nanmean(vort[wake_mask]**2))) if np.any(wake_mask) else 0.0
|
||||
|
||||
stable = nan_count == 0 and rho_min > 0.0 and rho_max < 2.0
|
||||
return {
|
||||
"nan_count": nan_count,
|
||||
"rho_min": rho_min,
|
||||
"rho_max": rho_max,
|
||||
"mass": mass,
|
||||
"ma_max": ma_max,
|
||||
"vort_rms": vort_rms,
|
||||
"stable": stable,
|
||||
}
|
||||
|
||||
|
||||
def plot_flow(rho, ux, uy, flag, nx, ny, title, out_path):
|
||||
"""Plot velocity magnitude, vorticity, and streamlines."""
|
||||
fluid_mask = flag.reshape(ny, nx) != FLUID
|
||||
vel = np.sqrt(ux**2 + uy**2)
|
||||
vel_m = np.ma.array(vel, mask=fluid_mask)
|
||||
vort = np.gradient(uy, axis=1) - np.gradient(ux, axis=0)
|
||||
vort_m = np.ma.array(vort, mask=fluid_mask)
|
||||
|
||||
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
|
||||
|
||||
# Velocity magnitude
|
||||
im0 = axes[0].imshow(vel_m, origin="lower", aspect="auto", cmap="turbo")
|
||||
plt.colorbar(im0, ax=axes[0], label="|u|")
|
||||
axes[0].set_title("Velocity Magnitude")
|
||||
|
||||
# Vorticity
|
||||
vals = vort[~fluid_mask]
|
||||
if vals.size > 0:
|
||||
vmax = max(float(np.percentile(np.abs(vals), 99)), 1e-8)
|
||||
else:
|
||||
vmax = 1e-6
|
||||
im1 = axes[1].imshow(vort_m, origin="lower", aspect="auto", cmap="RdBu_r",
|
||||
vmin=-vmax, vmax=vmax)
|
||||
plt.colorbar(im1, ax=axes[1], label="vorticity")
|
||||
axes[1].set_title("Vorticity")
|
||||
|
||||
# Streamlines
|
||||
X, Y = np.meshgrid(np.arange(nx), np.arange(ny))
|
||||
ux_s = np.ma.array(ux, mask=fluid_mask)
|
||||
uy_s = np.ma.array(uy, mask=fluid_mask)
|
||||
speed = np.ma.sqrt(ux_s**2 + uy_s**2)
|
||||
axes[2].streamplot(X, Y, ux_s, uy_s, color=speed, cmap="viridis",
|
||||
density=2.0, linewidth=0.7)
|
||||
axes[2].set_xlim(0, nx)
|
||||
axes[2].set_ylim(0, ny)
|
||||
axes[2].set_title("Streamlines")
|
||||
|
||||
fig.suptitle(title, fontsize=13)
|
||||
fig.tight_layout()
|
||||
fig.savefig(out_path, dpi=150)
|
||||
plt.close(fig)
|
||||
return out_path
|
||||
|
||||
|
||||
# ---------------------------------------------------------------------------
|
||||
# Case runner: double-buffer
|
||||
# ---------------------------------------------------------------------------
|
||||
def run_double_buffer(device_id, cfg, out_dir):
|
||||
"""Run a case with standard double-buffer streaming."""
|
||||
nx, ny = cfg["nx"], cfg["ny"]
|
||||
nq = cfg["nq"]
|
||||
n = nx * ny
|
||||
|
||||
set_macros(nx, ny, cfg["dim"], nq, cfg["vis"], cfg["u0"],
|
||||
cfg["collision_model"], cfg["use_les"], streaming_model=0,
|
||||
omega_collision_max=cfg.get("omega_max", 1.999),
|
||||
trt_magic_param=cfg.get("trt_magic", 0.1875))
|
||||
compiler.compile_kernel_v2()
|
||||
|
||||
cuda.init()
|
||||
dev = cuda.Device(device_id)
|
||||
ctx = dev.make_context()
|
||||
try:
|
||||
mod = cuda.module_from_file(compiler.kernel_path("kernel_v2.ptx"))
|
||||
init_fn = mod.get_function("InitTubeFlow_v2")
|
||||
step_fn = mod.get_function("OneStep")
|
||||
|
||||
# Upload d_params
|
||||
params_ptr, params_size = mod.get_global("d_params")
|
||||
params_data = pack_d_params(nx, ny, cfg["omega"], cfg["u0"])
|
||||
if len(params_data) < params_size:
|
||||
params_data += b"\x00" * (params_size - len(params_data))
|
||||
cuda.memcpy_htod(params_ptr, params_data)
|
||||
|
||||
fsize = n * nq * 4
|
||||
d_fi = cuda.mem_alloc(fsize)
|
||||
d_fi2 = cuda.mem_alloc(fsize)
|
||||
d_flag = cuda.mem_alloc(n)
|
||||
d_indx = cuda.mem_alloc(n * 4)
|
||||
d_delta = cuda.mem_alloc(4)
|
||||
d_action = cuda.mem_alloc(4)
|
||||
d_obs = cuda.mem_alloc(4)
|
||||
cuda.memset_d32(d_indx, 0, n)
|
||||
cuda.memset_d32(d_delta, 0, 1)
|
||||
cuda.memset_d32(d_action, 0, 1)
|
||||
cuda.memset_d32(d_obs, 0, 1)
|
||||
|
||||
block = (128, 1, 1)
|
||||
grid = ((nx + 127) // 128, ny, 1)
|
||||
|
||||
init_fn(d_flag, d_fi, block=block, grid=grid)
|
||||
cuda.memcpy_htod(d_flag, cfg["flag"])
|
||||
|
||||
host0 = np.empty(n * nq, dtype=np.float32)
|
||||
cuda.memcpy_dtoh(host0, d_fi)
|
||||
host0 = impose_rest_on_nonfluid(cfg["flag"], host0, nq, nx, ny)
|
||||
cuda.memcpy_htod(d_fi, host0)
|
||||
cuda.memcpy_htod(d_fi2, host0)
|
||||
|
||||
steps = cfg["steps"]
|
||||
report = max(steps // 5, 1)
|
||||
t0 = time.time()
|
||||
diverged_step = None
|
||||
for s in range(steps):
|
||||
step_fn(d_flag, d_fi, d_fi2, d_indx, d_delta, d_action, d_obs,
|
||||
block=block, grid=grid)
|
||||
d_fi, d_fi2 = d_fi2, d_fi
|
||||
|
||||
if (s + 1) % report == 0:
|
||||
cuda.Context.synchronize()
|
||||
h = np.empty(n * nq, dtype=np.float32)
|
||||
cuda.memcpy_dtoh(h, d_fi)
|
||||
rho_c = h.reshape(nq, ny, nx).sum(axis=0)
|
||||
nc = int(np.isnan(rho_c).sum())
|
||||
center = float(rho_c[ny // 2, nx // 2])
|
||||
print(f" step {s+1:6d}: rho_center={center:.6f} nan={nc}")
|
||||
if nc > 0:
|
||||
diverged_step = s + 1
|
||||
break
|
||||
|
||||
cuda.Context.synchronize()
|
||||
elapsed = time.time() - t0
|
||||
host = np.empty(n * nq, dtype=np.float32)
|
||||
cuda.memcpy_dtoh(host, d_fi)
|
||||
|
||||
rho, ux, uy = compute_macros_2d(host, nq, nx, ny, cfg["flag"])
|
||||
diag = diagnose(rho, ux, uy, cfg["flag"], nx, ny)
|
||||
diag["elapsed"] = elapsed
|
||||
diag["mlups"] = n * steps / elapsed / 1e6 if elapsed > 0 else 0
|
||||
diag["diverged_step"] = diverged_step
|
||||
|
||||
tag = cfg["tag"]
|
||||
plot_path = plot_flow(rho, ux, uy, cfg["flag"], nx, ny, tag,
|
||||
os.path.join(out_dir, f"{tag}.png"))
|
||||
diag["plot"] = plot_path
|
||||
return diag
|
||||
finally:
|
||||
ctx.pop()
|
||||
|
||||
|
||||
# ---------------------------------------------------------------------------
|
||||
# Case runner: Esoteric-Pull (single buffer)
|
||||
# ---------------------------------------------------------------------------
|
||||
def run_esopull(device_id, cfg, out_dir):
|
||||
"""Run a case with Esoteric-Pull single-buffer streaming."""
|
||||
nx, ny = cfg["nx"], cfg["ny"]
|
||||
nq = cfg["nq"]
|
||||
n = nx * ny
|
||||
|
||||
set_macros(nx, ny, cfg["dim"], nq, cfg["vis"], cfg["u0"],
|
||||
cfg["collision_model"], cfg["use_les"], streaming_model=1,
|
||||
omega_collision_max=cfg.get("omega_max", 1.999),
|
||||
trt_magic_param=cfg.get("trt_magic", 0.1875))
|
||||
compiler.compile_kernel_v2()
|
||||
|
||||
cuda.init()
|
||||
dev = cuda.Device(device_id)
|
||||
ctx = dev.make_context()
|
||||
try:
|
||||
mod = cuda.module_from_file(compiler.kernel_path("kernel_v2.ptx"))
|
||||
init_fn = mod.get_function("InitEsoPull")
|
||||
step_fn = mod.get_function("EsoPullStep")
|
||||
|
||||
# Upload d_params
|
||||
params_ptr, params_size = mod.get_global("d_params")
|
||||
params_data = pack_d_params(nx, ny, cfg["omega"], cfg["u0"])
|
||||
if len(params_data) < params_size:
|
||||
params_data += b"\x00" * (params_size - len(params_data))
|
||||
cuda.memcpy_htod(params_ptr, params_data)
|
||||
|
||||
fsize = n * nq * 4
|
||||
d_fi = cuda.mem_alloc(fsize)
|
||||
d_flag = cuda.mem_alloc(n)
|
||||
d_indx = cuda.mem_alloc(n * 4)
|
||||
d_delta = cuda.mem_alloc(4)
|
||||
d_action = cuda.mem_alloc(4)
|
||||
d_obs = cuda.mem_alloc(4)
|
||||
cuda.memset_d32(d_indx, 0, n)
|
||||
cuda.memset_d32(d_delta, 0, 1)
|
||||
cuda.memset_d32(d_action, 0, 1)
|
||||
cuda.memset_d32(d_obs, 0, 1)
|
||||
|
||||
block = (128, 1, 1)
|
||||
grid = ((nx + 127) // 128, ny, 1)
|
||||
|
||||
init_fn(d_flag, d_fi, block=block, grid=grid)
|
||||
cuda.memcpy_htod(d_flag, cfg["flag"])
|
||||
|
||||
# Note: for EsoPull, we don't impose_rest_on_nonfluid on the raw
|
||||
# DDF because the data is stored in esoteric layout. InitEsoPull
|
||||
# already stores rest equilibrium for solid nodes.
|
||||
|
||||
steps = cfg["steps"]
|
||||
report = max(steps // 5, 1)
|
||||
t0 = time.time()
|
||||
diverged_step = None
|
||||
|
||||
for s in range(steps):
|
||||
t_val = np.uint64(s) # timestep counter for load/store parity
|
||||
step_fn(d_fi, d_flag, d_indx, d_delta, d_action, d_obs,
|
||||
t_val, block=block, grid=grid)
|
||||
|
||||
if (s + 1) % report == 0:
|
||||
cuda.Context.synchronize()
|
||||
# For diagnostics, download raw DDF and decode from esopull layout
|
||||
h = np.empty(n * nq, dtype=np.float32)
|
||||
cuda.memcpy_dtoh(h, d_fi)
|
||||
# Esoteric layout: at this point the DDF is in post-store layout
|
||||
# for timestep s. To compute macros we need to "undo" the esoteric
|
||||
# read pattern. A simpler approach: compute rho = sum(fi) per node.
|
||||
# Because sum is invariant under slot permutation, rho is correct.
|
||||
# But ux/uy need correct direction assignment.
|
||||
# For diagnostic, use a simple sum-based stability check.
|
||||
f_arr = h.reshape(nq, ny, nx)
|
||||
rho_c = f_arr.sum(axis=0)
|
||||
nc = int(np.isnan(rho_c).sum())
|
||||
center = float(rho_c[ny // 2, nx // 2])
|
||||
print(f" step {s+1:6d}: rho_center={center:.6f} nan={nc}")
|
||||
if nc > 0:
|
||||
diverged_step = s + 1
|
||||
break
|
||||
|
||||
cuda.Context.synchronize()
|
||||
elapsed = time.time() - t0
|
||||
|
||||
# For final macros, do one more step that also writes to rho/u arrays.
|
||||
# But we don't have UpdateMacro for EsoPull yet. Instead, use the
|
||||
# approach: run a "read-only" macro computation from the esoteric layout.
|
||||
# For correctness, we load from the proper esoteric positions on host.
|
||||
h = np.empty(n * nq, dtype=np.float32)
|
||||
cuda.memcpy_dtoh(h, d_fi)
|
||||
|
||||
rho, ux, uy = _decode_esopull_macros(h, nq, nx, ny, cfg["flag"], steps)
|
||||
diag = diagnose(rho, ux, uy, cfg["flag"], nx, ny)
|
||||
diag["elapsed"] = elapsed
|
||||
diag["mlups"] = n * steps / elapsed / 1e6 if elapsed > 0 else 0
|
||||
diag["diverged_step"] = diverged_step
|
||||
|
||||
tag = cfg["tag"]
|
||||
plot_path = plot_flow(rho, ux, uy, cfg["flag"], nx, ny, tag,
|
||||
os.path.join(out_dir, f"{tag}.png"))
|
||||
diag["plot"] = plot_path
|
||||
return diag
|
||||
finally:
|
||||
ctx.pop()
|
||||
|
||||
|
||||
def _decode_esopull_macros(host_ddf, nq, nx, ny, flag, last_t):
|
||||
"""Decode macroscopic quantities from esoteric-pull layout on host.
|
||||
|
||||
After step t (0-based), the store was done at parity t.
|
||||
The next load would use parity t+1. To read correct DDFs we mimic
|
||||
load_f_esopull at t_read = last_t (the parity of the *next* step to execute).
|
||||
"""
|
||||
fi = host_ddf.reshape(nq, ny * nx) # fi[direction, node]
|
||||
t_read = last_t # parity for the load that would happen next
|
||||
|
||||
cx9 = np.array([0, 1, -1, 0, 0, 1, -1, 1, -1], dtype=np.float32)
|
||||
cy9 = np.array([0, 0, 0, 1, -1, 1, -1, -1, 1], dtype=np.float32)
|
||||
|
||||
# Compute neighbor table once
|
||||
j_table = np.zeros((nq, ny * nx), dtype=np.int64)
|
||||
for y in range(ny):
|
||||
for x in range(nx):
|
||||
k = y * nx + x
|
||||
xp = (x + 1) % nx
|
||||
xm = (x - 1) % nx
|
||||
yp = (y + 1) % ny
|
||||
ym = (y - 1) % ny
|
||||
j_table[0, k] = k
|
||||
j_table[1, k] = yp * nx + xp if nq > 1 else k # placeholder
|
||||
j_table[2, k] = ym * nx + xm if nq > 2 else k
|
||||
# D2Q9 neighbors: j[i] = neighbor in direction c_i
|
||||
if nq == 9:
|
||||
j_table[1, k] = y * nx + xp # +x
|
||||
j_table[2, k] = y * nx + xm # -x
|
||||
j_table[3, k] = yp * nx + x # +y
|
||||
j_table[4, k] = ym * nx + x # -y
|
||||
j_table[5, k] = yp * nx + xp # +x+y
|
||||
j_table[6, k] = ym * nx + xm # -x-y
|
||||
j_table[7, k] = ym * nx + xp # +x-y
|
||||
j_table[8, k] = yp * nx + xm # -x+y
|
||||
|
||||
n = nx * ny
|
||||
f_decoded = np.zeros((nq, n), dtype=np.float32)
|
||||
f_decoded[0] = fi[0]
|
||||
|
||||
for i in range(1, nq, 2):
|
||||
if t_read & 1:
|
||||
# Odd: f[i] from fi[n, i], f[i+1] from fi[j[i], i+1]
|
||||
f_decoded[i] = fi[i]
|
||||
f_decoded[i + 1] = fi[i + 1, j_table[i]]
|
||||
else:
|
||||
# Even: f[i] from fi[n, i+1], f[i+1] from fi[j[i], i]
|
||||
f_decoded[i] = fi[i + 1]
|
||||
f_decoded[i + 1] = fi[i, j_table[i]]
|
||||
|
||||
f_decoded = f_decoded.reshape(nq, ny, nx)
|
||||
rho = f_decoded.sum(axis=0)
|
||||
rho_safe = np.where(np.abs(rho) > 1e-12, rho, 1.0)
|
||||
ux = np.zeros_like(rho)
|
||||
uy = np.zeros_like(rho)
|
||||
for i in range(nq):
|
||||
ux += cx9[i] * f_decoded[i]
|
||||
uy += cy9[i] * f_decoded[i]
|
||||
ux /= rho_safe
|
||||
uy /= rho_safe
|
||||
return rho, ux, uy
|
||||
|
||||
|
||||
# ---------------------------------------------------------------------------
|
||||
# Case builders
|
||||
# ---------------------------------------------------------------------------
|
||||
def build_cases(steps_low, steps_high):
|
||||
"""Build the full test matrix."""
|
||||
# Grid params (moderate size for fast testing)
|
||||
nx, ny = 384, 192
|
||||
cx_ob, cy_ob, radius = 96.0, 96.0, 18.0
|
||||
u0 = 0.04
|
||||
|
||||
cases = []
|
||||
for re_val, re_label, n_steps, use_les in [
|
||||
(100.0, "Re100", steps_low, False),
|
||||
(100.0, "Re100", steps_low, True),
|
||||
(3000.0, "Re3000", steps_high, False),
|
||||
(3000.0, "Re3000", steps_high, True),
|
||||
]:
|
||||
for cm in (0, 1, 2):
|
||||
vis, omega = compute_vis_omega(re_val, 2.0 * radius, u0)
|
||||
les_tag = "LES" if use_les else "noLES"
|
||||
cm_name = COLLISION_NAMES[cm]
|
||||
tag = f"DB_{re_label}_{cm_name}_{les_tag}"
|
||||
cases.append({
|
||||
"tag": tag,
|
||||
"nx": nx, "ny": ny,
|
||||
"dim": 2, "nq": 9,
|
||||
"cx": cx_ob, "cy": cy_ob, "radius": radius,
|
||||
"flag": build_flags_2d(nx, ny, cx_ob, cy_ob, radius),
|
||||
"u0": u0,
|
||||
"vis": vis,
|
||||
"omega": omega,
|
||||
"collision_model": cm,
|
||||
"use_les": use_les,
|
||||
"steps": n_steps,
|
||||
"streaming": "double_buffer",
|
||||
"omega_max": 1.999,
|
||||
"trt_magic": 0.1875,
|
||||
})
|
||||
|
||||
# EsoPull case: low Re, SRT only
|
||||
re_eso = 100.0
|
||||
vis_eso, omega_eso = compute_vis_omega(re_eso, 2.0 * radius, u0)
|
||||
cases.append({
|
||||
"tag": "EsoPull_Re100_SRT_noLES",
|
||||
"nx": nx, "ny": ny,
|
||||
"dim": 2, "nq": 9,
|
||||
"cx": cx_ob, "cy": cy_ob, "radius": radius,
|
||||
"flag": build_flags_2d(nx, ny, cx_ob, cy_ob, radius),
|
||||
"u0": u0,
|
||||
"vis": vis_eso,
|
||||
"omega": omega_eso,
|
||||
"collision_model": 0,
|
||||
"use_les": False,
|
||||
"steps": steps_low,
|
||||
"streaming": "esopull",
|
||||
"omega_max": 1.999,
|
||||
"trt_magic": 0.1875,
|
||||
})
|
||||
|
||||
return cases
|
||||
|
||||
|
||||
# ---------------------------------------------------------------------------
|
||||
# Comparison plot: EsoPull vs DoubleBuffer
|
||||
# ---------------------------------------------------------------------------
|
||||
def plot_comparison(results, out_dir):
|
||||
"""Compare EsoPull and DoubleBuffer at matching Re/collision settings."""
|
||||
eso_key = "EsoPull_Re100_SRT_noLES"
|
||||
db_key = "DB_Re100_SRT_noLES"
|
||||
|
||||
eso = results.get(eso_key)
|
||||
db = results.get(db_key)
|
||||
if eso is None or db is None:
|
||||
return None
|
||||
|
||||
fig, axes = plt.subplots(2, 3, figsize=(18, 10))
|
||||
fig.suptitle("EsoPull vs DoubleBuffer — Re100 SRT noLES", fontsize=14)
|
||||
|
||||
labels = ["DoubleBuffer", "EsoPull"]
|
||||
for row, (r, label) in enumerate([(db, labels[0]), (eso, labels[1])]):
|
||||
vel_img = plt.imread(r["plot"]) if os.path.exists(r["plot"]) else None
|
||||
if vel_img is not None:
|
||||
axes[row, 0].imshow(vel_img)
|
||||
axes[row, 0].set_title(f"{label}: flow field")
|
||||
axes[row, 0].axis("off")
|
||||
else:
|
||||
axes[row, 0].text(0.5, 0.5, f"No image for {label}",
|
||||
ha="center", va="center", transform=axes[row, 0].transAxes)
|
||||
axes[row, 0].set_title(label)
|
||||
|
||||
# Metrics bar chart
|
||||
metrics = {
|
||||
"rho_min": r.get("rho_min", 0),
|
||||
"rho_max": r.get("rho_max", 0),
|
||||
"ma_max": r.get("ma_max", 0),
|
||||
"vort_rms": r.get("vort_rms", 0),
|
||||
}
|
||||
bars = list(metrics.keys())
|
||||
vals = [float(metrics[b]) for b in bars]
|
||||
axes[row, 1].barh(bars, vals, color=["steelblue", "salmon", "green", "purple"])
|
||||
axes[row, 1].set_title(f"{label}: diagnostics")
|
||||
|
||||
# Stability text
|
||||
text_lines = [
|
||||
f"stable: {r.get('stable', '?')}",
|
||||
f"nan_count: {r.get('nan_count', '?')}",
|
||||
f"mass: {r.get('mass', 0):.2f}",
|
||||
f"MLUPS: {r.get('mlups', 0):.1f}",
|
||||
f"diverged_step: {r.get('diverged_step', 'None')}",
|
||||
]
|
||||
axes[row, 2].text(0.1, 0.5, "\n".join(text_lines), fontsize=12,
|
||||
family="monospace", va="center",
|
||||
transform=axes[row, 2].transAxes)
|
||||
axes[row, 2].set_title(f"{label}: summary")
|
||||
axes[row, 2].axis("off")
|
||||
|
||||
fig.tight_layout()
|
||||
cmp_path = os.path.join(out_dir, "esopull_vs_doublebuffer.png")
|
||||
fig.savefig(cmp_path, dpi=150)
|
||||
plt.close(fig)
|
||||
return cmp_path
|
||||
|
||||
|
||||
# ---------------------------------------------------------------------------
|
||||
# Main
|
||||
# ---------------------------------------------------------------------------
|
||||
def main():
|
||||
parser = argparse.ArgumentParser(description="Stability matrix test")
|
||||
parser.add_argument("--device", type=int, default=0)
|
||||
parser.add_argument("--steps-low", type=int, default=3000,
|
||||
help="Steps for low-Re cases")
|
||||
parser.add_argument("--steps-high", type=int, default=6000,
|
||||
help="Steps for high-Re cases")
|
||||
parser.add_argument("--only-esopull", action="store_true",
|
||||
help="Only run the EsoPull test")
|
||||
args = parser.parse_args()
|
||||
|
||||
# Backup config/*.h files (kernel_v2.cu uses config.h, NOT macros.h)
|
||||
cfg_dir = os.path.join(os.path.dirname(compiler.kernel_path("config.h")), "config")
|
||||
config_files = ["config_grid.h", "config_physics.h", "config_method.h", "config_objects.h"]
|
||||
config_backups = {}
|
||||
for cf in config_files:
|
||||
path = os.path.join(cfg_dir, cf)
|
||||
with open(path, "r") as f:
|
||||
config_backups[path] = f.read()
|
||||
|
||||
out_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)),
|
||||
"..", "output", "stability_matrix")
|
||||
os.makedirs(out_dir, exist_ok=True)
|
||||
|
||||
cases = build_cases(args.steps_low, args.steps_high)
|
||||
if args.only_esopull:
|
||||
cases = [c for c in cases if c["streaming"] == "esopull"]
|
||||
|
||||
results = {}
|
||||
try:
|
||||
for i, cfg in enumerate(cases):
|
||||
tag = cfg["tag"]
|
||||
streaming = cfg["streaming"]
|
||||
print(f"\n[{i+1}/{len(cases)}] {tag}")
|
||||
print(f" Re={cfg['u0']*2*cfg['radius']/cfg['vis']:.0f}, "
|
||||
f"omega={cfg['omega']:.4f}, "
|
||||
f"collision={COLLISION_NAMES[cfg['collision_model']]}, "
|
||||
f"LES={cfg['use_les']}, streaming={streaming}")
|
||||
|
||||
if streaming == "esopull":
|
||||
diag = run_esopull(args.device, cfg, out_dir)
|
||||
else:
|
||||
diag = run_double_buffer(args.device, cfg, out_dir)
|
||||
|
||||
diag["tag"] = tag
|
||||
diag["streaming"] = streaming
|
||||
diag["collision"] = COLLISION_NAMES[cfg["collision_model"]]
|
||||
diag["use_les"] = cfg["use_les"]
|
||||
diag["re"] = cfg["u0"] * 2 * cfg["radius"] / cfg["vis"]
|
||||
results[tag] = diag
|
||||
|
||||
status = "PASS" if diag["stable"] else "FAIL"
|
||||
print(f" => {status}: rho=[{diag['rho_min']:.4f}, {diag['rho_max']:.4f}], "
|
||||
f"nan={diag['nan_count']}, ma_max={diag['ma_max']:.4f}, "
|
||||
f"MLUPS={diag['mlups']:.1f}")
|
||||
|
||||
# Comparison plot
|
||||
cmp_path = plot_comparison(results, out_dir)
|
||||
if cmp_path:
|
||||
print(f"\nComparison plot: {cmp_path}")
|
||||
|
||||
# Summary table
|
||||
print("\n" + "=" * 100)
|
||||
print(f"{'Tag':<35s} {'Stream':<8s} {'Col':<5s} {'LES':<5s} "
|
||||
f"{'Re':>6s} {'Stable':>7s} {'rho_min':>9s} {'rho_max':>9s} "
|
||||
f"{'Ma_max':>8s} {'MLUPS':>7s}")
|
||||
print("-" * 100)
|
||||
for tag, r in results.items():
|
||||
print(f"{tag:<35s} {r['streaming']:<8s} {r['collision']:<5s} "
|
||||
f"{'Y' if r['use_les'] else 'N':<5s} "
|
||||
f"{r['re']:6.0f} {'PASS' if r['stable'] else 'FAIL':>7s} "
|
||||
f"{r['rho_min']:9.5f} {r['rho_max']:9.5f} "
|
||||
f"{r['ma_max']:8.5f} {r['mlups']:7.1f}")
|
||||
print("=" * 100)
|
||||
|
||||
# Save JSON
|
||||
json_path = os.path.join(out_dir, "stability_matrix_results.json")
|
||||
json_results = {}
|
||||
for k, v in results.items():
|
||||
jr = {}
|
||||
for rk, rv in v.items():
|
||||
if isinstance(rv, (np.integer, np.floating)):
|
||||
jr[rk] = float(rv)
|
||||
elif isinstance(rv, np.bool_):
|
||||
jr[rk] = bool(rv)
|
||||
else:
|
||||
jr[rk] = rv
|
||||
json_results[k] = jr
|
||||
with open(json_path, "w") as f:
|
||||
json.dump(json_results, f, indent=2)
|
||||
print(f"\nResults saved: {json_path}")
|
||||
|
||||
finally:
|
||||
for path, content in config_backups.items():
|
||||
with open(path, "w") as f:
|
||||
f.write(content)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
Reference in New Issue
Block a user