泰勒-格林涡（收敛性研究）

Taylor-Green 涡是不可压缩流体的经典基准问题，具有已知的精确解：在正方形区域 \([0, 2\pi]^2\) 上的一个周期性衰减涡，其速度和压力均可写出闭式表达式。脚本 examples/fluid/taylor_green/taylor_green.py 利用该问题对瞬态 Taylor-Hood Navier-Stokes 求解器进行定量的 h-收敛性研究 —— 这是 流体力学 中唯一以验证而非可视化为目的的算例。

精确解

对于黏度 \(\nu\) 与时间 \(t\)，

\[\begin{split}u(x, y, t) &= -\cos(x)\, \sin(y)\, e^{-2\nu t}, \\ v(x, y, t) &= \phantom{-}\sin(x)\, \cos(y)\, e^{-2\nu t}, \\ p(x, y, t) &= -\tfrac14 \bigl(\cos(2x) + \cos(2y)\bigr)\, e^{-4\nu t}.\end{split}\]

速度以速率 \(2\nu\)（运动黏度的两倍）指数衰减；压力以速率 \(4\nu\) 衰减。代入即可验证 \((u, v, p)\) 精确满足不可压缩纳维-斯托克斯方程。

该脚本在每个边界节点上施加等于当前时刻精确解的狄利克雷速度，并将一个压力自由度固定为其精确值。这避免了周期性边界条件的需求，并分离出离散化误差。

求解器

直接复用来自 圆柱绕流（涡脱落） 的瞬态 NavierStokesTransientAssembler —— 空间上采用 Taylor-Hood P2-P1（MixedElementAssembler），时间上采用后向欧拉加 Picard 子迭代，forward_vector 用于 \(\rho/\Delta t\,\mathbf{u}^n\cdot\mathbf{v}\) 载荷。唯一的改动如下：

  • 计算域为 \([0, 2\pi]^2\)（二阶生成网格，因此 P2 速度节点即为网格点，精确解可直接在其上求值），

  • 狄利克雷速度取自精确解而非抛物线入口 —— 每步通过 condenser.update_dirichlet 更新，

  • 最后时刻速度及压力的 \(L^2\) 误差与解析解对比。

列表 24 examples/fluid/taylor_green/taylor_green.py（核心部分）
pin = layout.dof_index("p", int(layout.node_ids("p")[0]))

def dirichlet_values(t):
    bc_val = torch.zeros(layout.n_dofs, dtype=torch.float64)
    bc_val[layout.dof_mask("u")] = exact_velocity(points, t, nu).reshape(-1)
    bc_val[pin] = exact_pressure(x_p, t, nu)[0]
    return bc_val

for step in range(1, n_steps + 1):
    condenser.update_dirichlet(dirichlet_values(step * dt))
    u_prev = layout.split(sol)["u"]
    for _ in range(picard_iter):
        w = layout.split(u_iter)["u"]
        K = assembler(point_data={"w": w})
        f = assembler.assemble_vector(point_data={"uprev": u_prev})
        K_, f_ = condenser(K, f)
        u_iter = condenser.recover(K_.solve(f_))
    sol = u_iter

h 收敛性

Taylor-Hood 空间理论预示 \(\|\mathbf{u}_h - \mathbf{u}\|_{L^2} = \mathcal{O}(h^3)\) 且 \(\|p_h - p\|_{L^2} = \mathcal{O}(h^2)\) —— 但后向欧拉额外引入 \(\mathcal{O}(\Delta t)\) 误差。因此该研究取 \(\Delta t = h^2/4\)，使得时间误差至少以与 \(\mathcal{O}(h^2)\) 压力误差同等的速度衰减。在此缓慢衰减的涡流（\(\nu = 0.01\)）中，时间误差常数较小，故实测收敛率落在 空间 极限 —— 速度约为 \(3\)，P1 压力约为 \(2\)：

Grid   h          L2_vel         L2_pres        Rate_vel   Rate_pres
----------------------------------------------------------------
10     0.6283     1.423434e-01   3.209651e-01   -          -
20     0.3142     1.862684e-02   8.626960e-02   2.93       1.90
40     0.1571     1.434308e-03   2.045441e-02   3.70       2.08

taylor_green_convergence.png 以双对数坐标展示相同数据，并标注参考 \(h^2\) 斜率。若观测收敛率降至 1，则求解器中存在错误 —— 最常见的是边界条件施加或时间步进中的隐蔽缺陷。

质量加权误差范数

该脚本通过 P2 空间的质量矩阵计算离散 \(L^2\) 范数：

\[\|e_h\|_{L^2}^2 \;=\; e_h^T\, M\, e_h \;=\; \sum_K \int_K e_h^2 \,\mathrm{d}\Omega,\]

由 MassElementAssembler 在二阶网格上装配。P1 压力先延拓至 P2 空间（layout.prolong —— 精确的，因为 \(P_1 \subset P_2\)），故一个质量矩阵即可服务于两个场。这是有限元误差分析中的正确范数，不同于节点误差向量的简单欧几里得范数。

输出

  • 控制台表格——每个加密步上的收敛阶。

  • ``taylor_green_convergence.png`` —— 误差双对数图随网格尺寸变化，附参考斜率。

  • ``taylor_green_results.png``——最细网格的三联面板快照：速度、压力、速度误差大小。

  • ``taylor_green.mp4``（可选）——衰减涡的动画。

t=0.5 时的泰勒-格林涡——涡量＋流线，速度＋速度矢量

图 65 taylor_green.py 在 \(t=0.5\) 时的输出。左：叠加流线的涡量场——经典的周期性 \(2\times2\) 交替符号涡阵列。右：带速度矢量的速度大小，展示了相邻卷之间特征性的“鞍”状结构。相对解析初始条件，幅值已衰减了 \(\exp(-2\nu t)\)，脚本的收敛性研究即以此为参考真值。

运行方式

cd examples/fluid/taylor_green
python taylor_green.py      # writes convergence + results pngs

下一步