圆柱绕流（涡脱落）

稳态方腔算例的瞬态对应版本。脚本 examples/fluid/cylinder_flow/cylinder_flow.py 运行经典的圆柱绕流基准——一个细长的矩形通道，入口附近放置一个小圆柱。在 \(\mathrm{Re} = 100\) 时尾迹变得不稳定，涡从圆柱的上下两侧交替脱落，一条冯·卡门涡街向下游传播。其几何与参数遵循 DFG 2D-1 基准（Schäfer & Turek，1996），该基准给出了阻力/升力系数与斯特劳哈尔数的参考值。

问题

瞬态不可压缩纳维-斯托克斯方程组，

\[\rho\, \frac{\partial \mathbf{u}}{\partial t} + \rho\, (\mathbf{u} \cdot \nabla)\mathbf{u} \;=\; -\nabla p + \mu\, \Delta \mathbf{u}, \qquad \nabla \cdot \mathbf{u} = 0,\]

定义在通道 \(\Omega = [0, 2.2] \times [0, 0.41]\) 上，其中含有一个半径 \(r = 0.05\)、圆心位于 \((0.2, 0.2)\) 的圆柱。边界条件：

  • 入口（\(x = 0\)）：抛物型分布 \(u_x(y) = 4\,U_\text{max}\, y\, (H - y) / H^2\)，\(u_y = 0\)，

  • 壁面与圆柱表面：无滑移，

  • 出口（\(x = 2.2\)）：“自由出流”（do-nothing）＋压力钉定，

取 \(U_\text{max} = 1.5\) 得到 \(\bar{U} = 1\)，\(D = 0.1\)，\(\rho = 1\)，\(\mu = 10^{-3}\)，由此 \(\mathrm{Re} = \rho \bar{U} D / \mu = 100\)。

时间积分：隐式欧拉法

瞬态项 \(\rho\, \partial_t \mathbf{u}\) 采用后向欧拉法离散——隐式且无条件稳定。每个时间步求解

\[\rho\, \frac{\mathbf{u}^{n+1} - \mathbf{u}^{n}}{\Delta t} + \rho\, (\mathbf{u}^{n} \cdot \nabla)\mathbf{u}^{n+1} \;=\; -\nabla p^{n+1} + \mu\, \Delta \mathbf{u}^{n+1},\]

（对流项通过将输运速度滞后取为 \(\mathbf{u}^n\) 来线性化）。脚本将其实现为方腔装配器的瞬态变体，在对角块上加入质量项 \(\rho/\Delta t \cdot u\, v\) 以及相应的 SUPG 残差：

列表 16 examples/fluid/cylinder_flow/cylinder_flow.py（核心部分）
class NavierStokesTransientAssembler(ElementAssembler):
    def __post_init__(self, rho=1.0, mu=0.01, dt=1e-3, tau=1e-3):
        self.rho, self.mu, self.dt, self.tau = rho, mu, dt, tau

    def forward(self, u, v, gradu, gradv, w_prev):
        dim = gradu.shape[0]
        mass        = self.rho * (u * v) / self.dt
        convection  = self.rho * torch.dot(w_prev, gradv) * u
        diffusion   = self.mu * torch.dot(gradu, gradv)
        supg_test   = self.tau * torch.dot(w_prev, gradu)
        supg_res    = (self.rho/self.dt * v
                       + self.rho * torch.dot(w_prev, gradv)) * supg_test
        momentum_diag = mass + convection + diffusion + supg_res
        # …assemble (dim+1)×(dim+1) block as in cavity.py…

上一时间步的解 \(\mathbf{u}^n\) 通过一个独立的 NodeAssembler（MomentumRHSAssembler）进入右端项，贡献 \(\rho/\Delta t\, \mathbf{u}^n \cdot v\)。两个装配器共享同一个 tau 场，每一步都会更新以反映局部网格单元的雷诺数。

稳定化

圆柱问题比方腔更为苛刻：上游网格单元看到的几乎是均匀流，而尾迹区域则需要精细的对流主导型稳定化。脚本采用一个自适应的 \(\tau\)：

\[\tau \;=\; \left[ \left(\frac{2}{\Delta t}\right)^2 + \left(\frac{2|\mathbf{u}|}{h}\right)^2 + \left(\frac{4\nu}{h^2}\right)^2 \right]^{-1/2},\]

它在每个单元上由局部网格尺寸 \(h\) 与速度大小计算得到。这使得 SUPG 的扩散部分在流速较慢处占主导，而对流部分在尾迹中占主导。

后处理

在每个保存的帧上，脚本通过 \(L^2\) 投影导出一个涡量场 \(\omega = \partial_x v - \partial_y u\)：装配一个质量矩阵 MassElementAssembler，由 \(\partial_x v - \partial_y u\)（其本身用 NodeAssembler 装配）构造投影右端项 \(\tilde\omega\)，再求解 \(M\,\omega = \tilde\omega\)。这是在分片线性有限元上从逐单元梯度恢复光滑节点场的标准做法。

冯·卡门涡街是需要关注的定性特征：一旦尾迹失稳（通常在 \(t \gtrsim 5\) 之后），涡便以斯特劳哈尔数 \(\mathrm{St} = f D / \bar{U} \approx 0.3\) 从圆柱的上下两侧交替脱落，可在涡量动画中直接观察到。

输出与渲染

  • 帧序列。每隔 N 步，脚本通过 mesh.plot 渲染一张三联面板 PNG（涡量、速度、压力）到 frames/ 中。

  • MP4 渲染。配套脚本 examples/fluid/cylinder_flow/render_video.py 通过一次 ffmpeg 拼接将 frames/*.png 序列合成为 vortex_street.mp4。

  • 最终快照。cylinder_flow_final.png 是最后一步的同款三联面板图，便于用于演讲与报告。

``cylinder_flow.py`` 的输出（由 render_video.py 渲染为 MP4）：圆柱后方发展完整的冯·卡门涡街。在初始的对称阶段之后，一个微小的非对称扰动触发周期性脱落；涡的符号交替变化，并以大致等于入口速度的速度向下游输运。

运行方式

cd examples/fluid/cylinder_flow
python cylinder_flow.py            # writes frames/*.png + cylinder_flow_final.png
python render_video.py             # stitches frames/ into vortex_street.mp4

这一瞬态运行是整个示例库中耗时最长的——默认配置为数千个时间步。若想快速冒烟测试，可减小 n_steps 或粗化网格。

下一步