圆柱绕流（涡脱落）

稳态空腔算例的瞬态对应版本。脚本 examples/fluid/cylinder_flow/cylinder_flow.py 运行经典的圆柱绕流基准算例——一个长矩形通道，在入口附近放置一个小圆柱。当 \(\mathrm{Re} = 100\) 时，尾流失稳，涡旋交替从圆柱上下表面脱落，卡门涡街向下游传播。几何与参数遵循 DFG 二维基准（Schäfer & Turek, 1996）。

离散化：通过 MixedElementAssembler 采用 Taylor-Hood P2-P1 元，时间方向用向后欧拉格式，配合 Picard 子迭代。以下两点使该脚本成为"进阶"混合元示例：

  • gmsh/MeshGen 生成的通道网格是线性的——二次速度空间在其上拓扑地构造（每条唯一边对应一个自由度），无需二阶重新剖分；

  • 载荷向量与涡量后处理均通过 assemble_vector() 完成，滞后速度借助 field_data 搭载在 P2 场自身的自由度上。

问题

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

\[\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\)。

时间积分：向后欧拉 + Picard

瞬态项采用向后欧拉离散——隐式、无条件稳定。每个时间步求解

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

其中对流速度 \(\mathbf{w}\) 每步通过一到两次 Picard 子迭代更新。矩阵端与载荷端由同一个装配器的两部分构成：

列表 21 examples/fluid/cylinder_flow/cylinder_flow.py（核心部分）
class NavierStokesTransientAssembler(MixedElementAssembler):
    fields = [
        Field(trial="u", test="v", order=2, components=2),
        Field(trial="p", test="q", order=1),
    ]

    def __post_init__(self, rho=1.0, mu=0.01, dt=1e-3):
        self.rho, self.mu, self.dt = rho, mu, dt

    def forward(self, u, gradu, p, v, gradv, q, w):
        mass = self.rho / self.dt * u.dot(v)
        convection = self.rho * (gradu @ w).dot(v)
        diffusion = self.mu * (gradu * gradv).sum()
        return mass + convection + diffusion \
            - p * gradv.diagonal().sum() \
            - q * gradu.diagonal().sum()

    def forward_vector(self, v, uprev):
        return self.rho / self.dt * uprev.dot(v)

时间循环中每次 Picard 迭代仅需三行装配代码——注意 field_data 通道承载滞后的 P2 速度：

for step in range(n_steps):
    u_prev = layout.split(sol)["u"]        # [n_u, 2] on the P2 DOFs
    for _ in range(picard_iter):
        w = layout.split(u_iter)["u"]
        K = assembler(field_data={"w": ("u", w)})
        f = assembler.assemble_vector(field_data={"uprev": ("u", u_prev)})
        K_, f_ = condenser(K, f)
        u_iter = condenser.recover(K_.solve(f_))
    sol = u_iter

拓扑 P2 空间上的边界条件

MeshGen 网格不携带 is_boundary 点数据，且半数速度节点是根本不是网格点的边中点。两个问题均可通过布局的拓扑辅助函数解决——boundary_mask 按面关联关系对自由度分类，points("u") 为每个速度节点赋予坐标：

x_u = layout.points("u")
is_boundary = layout.split(layout.boundary_mask("u"))["u"][:, 0]
is_inlet = is_boundary & (x_u[:, 0] <= eps)
is_outlet = x_u[:, 0] >= length - eps
no_slip = is_boundary & ~is_inlet & ~is_outlet   # walls + cylinder

bc_mask = layout.dof_mask("u", node_mask=is_inlet | no_slip)
y_in = x_u[is_inlet, 1]
bc_val[layout.dof_mask("u", node_mask=is_inlet, component=0)] = \
    4.0 * u_max * y_in * (height - y_in) / (height * height)

出口保持自由（do-nothing），从而锚定压力——鞍点系统无需点约束即可非奇异。

后处理：以混合载荷向量形式提取涡量

每个保存帧时刻，脚本通过 \(L^2\) 投影将涡量 \(\omega = \partial_x u_y - \partial_y u_x\) 恢复到 P1 压力空间。投影右端项 \(\int \omega_h\, q\,\mathrm{d}x\) 采用速度的精确 P2 梯度——一行 func= 线性形式即可，速度通过 field_data 传入：

omega_rhs = assembler.assemble_vector(
    func=lambda q, gradw: (gradw[1, 0] - gradw[0, 1]) * q,
    field_data={"w": ("u", velocity)},
)
omega = m_mat.solve(layout.split(omega_rhs)["p"])   # P1 mass matrix

卡门涡街是应关注的定性特征：一旦尾流失稳，涡旋即交替从圆柱上下表面脱落，并以近似平均入口速度向下游对流。

输出与渲染

  • 帧序列。 每隔 save_every 步，脚本通过 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 或粗化网格。

下一步