顶盖驱动方腔¶

顶盖驱动方腔流是不可压缩 Navier-Stokes 求解器的经典基准问题。其物理机制十分简单——一个充满流体的方腔，顶壁以单位速度滑动，其余壁面为无滑移边界，无体积力——但在中等雷诺数下，流场已呈现主涡及次级角涡结构，任何合理的求解器都必须能够复现这些特征。examples/fluid/cavity/ 目录下的两个脚本采用 Taylor-Hood P2-P1 混合单元求解定常问题，其中 cavity.py 在三角形剖分的单位方域上运行，cavity_3d.py 则在四面体剖分的单位立方体上运行。二者共享相同的标量弱形式被积函数——该被积函数与空间维度无关，因此三维情形本质上只需替换网格并设置 components=3。

由于 P2-P1 单元对满足离散 inf-sup（LBB）条件，标准 Galerkin 形式本身即稳定：这些脚本中不存在任何 SUPG/PSPG 稳定化项，也不含 tau 参数。块 bookkeeping（两个场、不同阶次、不同节点集）由 MixedElementAssembler 负责管理——有关机制详见混合装配，线性（Stokes）预热算例参见 Taylor-Hood Stokes。

问题¶

稳态不可压缩纳维-斯托克斯方程的强形式为

\[\rho\, (\mathbf{u} \cdot \nabla)\mathbf{u} \;=\; -\nabla p + \mu\, \Delta \mathbf{u} \quad \text{in } \Omega, \qquad \nabla \cdot \mathbf{u} \;=\; 0,\]

其中 \(\Omega = (0, 1)^d\)、\(\mu = 1/\mathrm{Re}\)、\(\rho = 1\)，边界条件为

顶盖（\(y = 1\)）：\(\mathbf{u} = (1, 0, \dots)\)（运动壁面）
其余壁面：\(\mathbf{u} = \mathbf{0}\)（无滑移）
一个压力自由度被固定为零（方腔为封闭区域，压力仅确定到相差一个常数）。

弱形式与 Picard 线性化¶

以试函数 \((\mathbf{v}, q)\) 相乘，并将对流速度滞后至前一次 Picard 迭代 \(\mathbf{w} = \mathbf{u}^{n}\)，得到线性化混合形式：求 \((\mathbf{u}, p)\) 使得

\[\int_\Omega \rho\,(\mathbf{w}\cdot\nabla)\mathbf{u}\cdot\mathbf{v}\,\mathrm{d}x + \int_\Omega \mu\,\nabla\mathbf{u} : \nabla\mathbf{v}\,\mathrm{d}x - \int_\Omega p\,(\nabla\cdot\mathbf{v})\,\mathrm{d}x - \int_\Omega q\,(\nabla\cdot\mathbf{u})\,\mathrm{d}x \;=\; 0\]

对所有 \((\mathbf{v}, q)\) 成立。采用 Taylor-Hood 空间无需额外处理，TensorMesh 的实现是直接转录——两个场声明及标量被积函数如下：

列表 19 examples/fluid/cavity/cavity.py（核心部分）¶

class NavierStokesAssembler(MixedElementAssembler):
    fields = [
        Field(trial="u", test="v", order=2, components=2),  # P2 velocity
        Field(trial="p", test="q", order=1),                # P1 pressure
    ]

    def __post_init__(self, rho=1.0, mu=0.01):
        self.rho = rho
        self.mu = mu

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

gradu 为速度雅可比矩阵（二维时为 [2, 2]，三维时为 [3, 3]），因此 (gradu @ w).dot(v) 对应 convection 项 \((\mathbf{w}\cdot\nabla)\mathbf{u}\cdot\mathbf{v}\)，而 gradu.diagonal().sum() 为散度。滞后速度 w 为数据（非试/检验场），故被积函数保持双线性——这正是混合装配器的约定。由于表达式未硬编码维度，同一类即可服务于二维与三维。

Picard 迭代¶

每次迭代以当前速度迭代值重新装配矩阵，并求解线性鞍点系统：

mesh = Mesh.gen_rectangle(chara_length=1.0 / n_grid, order=2).double()
assembler = NavierStokesAssembler.from_mesh(mesh, rho=1.0, mu=1.0 / re)
layout = assembler.layout

bc_mask = layout.dof_mask("u", mesh.boundary_mask)   # no-slip everywhere
bc_mask[layout.dof_index("p", int(layout.node_ids("p")[0]))] = True  # pin

bc_val = torch.zeros(layout.n_dofs, dtype=torch.float64)
bc_val[layout.dof_mask("u", is_top, component=0)] = 1.0   # moving lid

condenser = Condenser(bc_mask, bc_val[bc_mask])
sol = torch.zeros(layout.n_dofs, dtype=torch.float64)
sol[bc_mask] = bc_val[bc_mask]

for i in range(max_iter):
    w = layout.split(sol)["u"]            # previous-iterate velocity
    K = assembler(point_data={"w": w})
    f = torch.zeros(layout.n_dofs, dtype=torch.float64)

    K_, f_ = condenser(K, f)
    sol = condenser.recover(K_.solve(f_))
    # ...convergence check on the relative update...

几个需要注意的细节：

块自由度布局。 解向量先将所有速度自由度堆叠在前，随后是压力自由度。layout.dof_mask / layout.dof_index 按该编号构建边界掩码，layout.split(sol) 返回 {"u": [n_u, 2], "p": [n_p]}——无需任何手写的索引运算。
顶盖掩码。 在二阶网格上，P2 速度节点恰好与网格点重合，因此 point_data={"w": w} 及节点掩码如 is_top = mesh.points[:, 1] > 1 - 1e-6 可直接应用。
求积。 默认阶次（2 * max(order) = 4）对 convection 项的积分精度差一阶——这是标准且无害的选择；如需精确积分，可传入 quadrature_order=。

Picard 迭代在中等雷诺数下呈几何收敛——在默认 \(30 \times 30\) 网格上，Re = 100 时相对更新量在 8 次迭代内降至 \(10^{-4}\) 以下。

Re=100 时顶盖驱动方腔的速度大小与压力 — 图 61 `cavity.py` 在 Re = 100（Taylor-Hood P2-P1）下的输出。左：速度模长 \(\|u\|\)——运动顶盖将流体拖向右上方，沿右壁下行并形成主涡。右：P1 压力（延拓至 P2 网格点以绘制），右上角存在特征性高压区，该处顶盖与壁面滞止。¶

扩展到 3D——`cavity_3d.py`¶

cavity_3d.py 在单位立方体上求解相同物理问题，展示了获取 Taylor-Hood 单元对的第二种方式：Mesh.gen_cube 生成的四面体网格为线性的，二次速度空间由混合装配器 拓扑地 生成（四面体网格每条唯一边附加一个额外自由度）——无需二阶重网格。场声明变为

fields = [
    Field(trial="u", test="v", order=2, components=3),
    Field(trial="p", test="q", order=1),
]

forward 被积函数与二维情形逐字节相同。由于速度自由度不再是网格点，脚本使用布局的拓扑辅助函数处理边界条件与后处理：

列表 20 examples/fluid/cavity/cavity_3d.py（核心部分）¶

x_u = layout.points("u")                                   # P2 node coordinates
is_boundary = layout.split(layout.boundary_mask("u"))["u"][:, 0]
is_top = is_boundary & (x_u[:, 1] > 1.0 - 1e-6)

bc_mask = layout.dof_mask("u", node_mask=is_boundary)      # no-slip walls
bc_val[layout.dof_mask("u", node_mask=is_top, component=0)] = 1.0

# Picard loop: the lagged velocity lives on the P2 field's own DOFs
K = assembler(field_data={"w": ("u", w)})

# post-processing: interpolate the P2 velocity back to mesh points
velocity = layout.prolong("u", layout.split(sol)["u"])

注意与二维的两处差异：boundary_mask 以拓扑方式寻找壁面（面关联——边生成节点无需 is_boundary 点数据），且滞后速度通过 field_data 传入，因其无网格点表示。在 chara_length=0.1 时，系统含 24k 自由度（1.1k 网格点上的 7.6k P2 速度节点），13 次 Picard 迭代收敛。输出为体数据：cavity_3d.vtu 供 ParaView 使用，另含通过 PyVista 渲染的 \(z = 0.5\) 中截面切片。

3D 顶盖驱动方腔在 z=0.5 中平面上的速度与压力 — 图 62 `cavity_3d.py` 输出：立方体 \(z=0.5\) 中截面切片上的速度模长（左）与压力（右）。流型在顶盖附近与二维解吻合，但向前、后壁衰减，故中截面显示较二维基准更弱的回流。¶

运行方式¶

cd examples/fluid/cavity
python cavity.py        # 2D, writes cavity_results.png
python cavity_3d.py     # 3D, writes cavity_3d.vtu + cavity_3d.png

在默认分辨率下，二者于 Re = 100 时均在 \(\lesssim 15\) 次 Picard 迭代内收敛。

下一步¶

Taylor-Hood Stokes — 含收敛性研究的线性预热算例。
圆柱绕流（涡脱落） — 瞬态版本：向后 Euler、forward_vector 载荷、涡脱落。
瑞利-贝纳德对流 — 向同一块系统添加第三个场（温度）。
混合装配 — 混合装配约定的完整说明。