瑞利-贝纳德对流

瑞利-贝纳德问题将不可压缩流体流动与热输运耦合：一个水平空腔从底部加热、顶部冷却，当浮力超过粘性和热扩散的耗散效应时，会形成对流卷。脚本 examples/fluid/rayleigh_benard/rayleigh_benard.py 在宽高比为 2:1 的矩形空腔中求解布辛涅斯克近似，瑞利数可配置。

这是图库中的 三场混合 示例：泰勒-胡德 P2-P1 速度/压力 外加 P2 温度，全部声明在一个 MixedElementAssembler 上。浮力耦合和能量方程只是同一标量被积函数的更多项——无需逐节点块填充，无需稳定化。

问题

瞬态布辛涅斯克方程：

\[ \begin{align}\begin{aligned}\rho\,\frac{\partial\mathbf{u}}{\partial t} + \rho\, (\mathbf{u} \cdot \nabla)\mathbf{u} \;=\; -\nabla p + \mu\, \Delta \mathbf{u} + \rho\, g\, \beta\, T\, \hat{\mathbf{e}}_y, \qquad \nabla \cdot \mathbf{u} = 0,\\\frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla T \;=\; \kappa\, \Delta T,\end{aligned}\end{align} \]

定义在 \(\Omega = [0, 2] \times [0, 1]\) 上，其中

  • 速度：所有壁面均为无滑移，

  • 温度：底部 (\(y = 0\)) \(T = 1\)，顶部 (\(y = 1\)) \(T = 0\)，侧壁无热通量，

  • 一个压力自由度固定（封闭空腔）。

瑞利数

\[\mathrm{Ra} \;=\; \frac{g\, \beta\, \Delta T\, L^3}{\nu\, \alpha}\]

控制流动状态：当瑞利数低于 \(\mathrm{Ra}_c \approx 1708\) 时，传热为纯导热；超过该值则出现对流卷。脚本默认使用 \(\mathrm{Ra} = 2 \times 10^4\)（\(\mathrm{Pr} = 1\)），处于稳定的对流状态。

为何进行时间推进？ 导热状态（线性温度分布、零速度）在 任意 瑞利数下都是稳态解——从静止开始的稳态求解器会直接收敛到该解。当瑞利数超过 \(\mathrm{Ra}_c\) 时，该状态不稳定，因此脚本对导热温度分布施加微小扰动，并采用向后欧拉法进行时间积分：不稳定性会物理地增长为对流卷，当解不再变化时运行停止。

三场，一个被积函数

未知量声明为三个场——试函数名 (u, p, T)，测试函数名 (v, q, s)。耦合弱形式的每一项将一个试函数因子与一个测试函数因子配对，因此整个系统是一个双线性被积函数：

列表 23 examples/fluid/rayleigh_benard/rayleigh_benard.py（核心部分）
class RayleighBenardAssembler(MixedElementAssembler):
    fields = [
        Field(trial="u", test="v", order=2, components=2),  # P2 velocity
        Field(trial="p", test="q", order=1),                # P1 pressure
        Field(trial="T", test="s", order=2),                # P2 temperature
    ]

    def __post_init__(self, rho=1.0, mu=0.1, kappa=0.1,
                      g=10.0, beta=1.0, dt=1e-2):
        self.rho, self.mu, self.kappa = rho, mu, kappa
        self.g, self.beta, self.dt = g, beta, dt

    def forward(self, u, gradu, p, T, gradT, v, gradv, q, s, grads, w):
        momentum = self.rho / self.dt * u.dot(v) \
            + self.rho * (gradu @ w).dot(v) \
            + self.mu * (gradu * gradv).sum() \
            - p * gradv.diagonal().sum() \
            - self.rho * self.g * self.beta * T * v[1]   # buoyancy
        continuity = -q * gradu.diagonal().sum()
        energy = T * s / self.dt \
            + w.dot(gradT) * s \
            + self.kappa * gradT.dot(grads)
        return momentum + continuity + energy

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

两种耦合值得仔细分析：

  • 浮力项 - rho g beta T v[1] 将 试函数温度 与 垂直速度测试函数 配对——一个非对角 \((v, T)\) 块，由混合装配器自动提取。因此温度以 全隐式 方式反馈到动量方程中；只有对流速度 w 采用滞后处理（皮卡迭代）。

  • 能量输运 w.dot(gradT) * s + kappa * gradT.dot(grads) 是嵌在同一矩阵中的标量对流-扩散方程，采用相同的滞后速度 w。

时间推进至吸引子

每个向后欧拉步通过同一装配器，用当前速度装配矩阵，用上一时刻的场装配载荷向量：

T0 = (1.0 - y) + 0.01 * torch.sin(math.pi * x / 2) * torch.sin(math.pi * y)
sol = layout.cat(u=0.0, p=0.0, T=T0)

for step in range(n_steps):
    fields = layout.split(sol)
    K = assembler(point_data={"w": fields["u"]})
    f = assembler.assemble_vector(
        point_data={"uprev": fields["u"], "Tprev": fields["T"]})
    K_, f_ = condenser(K, f)
    sol = condenser.recover(K_.solve(f_))
    # stop when the per-step relative update stalls (steady state)

在二阶网格上，两个 P2 场（速度和温度）均定义在网格节点上，因此滞后数据通过普通 point_data 传递，最终场可直接绘图。

瑞利-贝纳德对流的温度与速度大小

图 64 rayleigh_benard.py 的输出。左：温度场——温暖的底部壁面产生上升羽流，在冷顶部分裂为向下流动的冷羽流。右：速度大小——反向旋转的对流卷，上升/下沉柱处速度最大，卷心处为滞止点。

运行方式

cd examples/fluid/rayleigh_benard
python rayleigh_benard.py     # writes rayleigh_benard.png

编辑 ra= 参数以扫描瑞利数；远高于 \(10^5\) 的值会变得非稳态，需要更细的网格和更小的时间步。

下一步