网格一歪，拉普拉斯算子就漏 — 非正交扩散通量修正

立方体上跑得好好的扩散求解器,搬到歪斜的三角形上,答案会悄悄偏掉。方程相同。代码也相同。温度场变糊,收敛变慢。不是边界条件,也不是物性。是面上量梯度的方式。

把非正交角滑块从20°推到65°,过松弛的青色矢量会长过面积,最小修正的粉色矢量拉得更长。

青色比粉色短的时候,画面下方的 ρ 越过 1。看45°时是不是只有最小修正先变红。

面上的梯度才是麻烦

有限体积法(FVM,按单元积分守恒量)里,扩散项是穿过面的通量之和。面 $f$ 上的扩散通量这样写。

$\Gamma_f$ 是面上的扩散系数。$\nabla\phi_f$ 是面梯度。$\mathbf{S}_f$ 是外向面积矢量,面积乘法向。

难的是把 $\nabla\phi_f \cdot \mathbf{S}_f$ 收成单元值。正交网格上,连接两个单元中心的 $\mathbf{d}$ 与法向平行。$(\phi_N-\phi_P)/|\mathbf{d}|$ 就是法向梯度。

非结构网格上这个假设破了。$\mathbf{d}$ 和 $\mathbf{S}_f$ 错开。那个角是非正交角 $\theta$。$(\phi_N-\phi_P)/|\mathbf{d}|$ 现在是沿 $\mathbf{d}$ 的梯度,不是法向梯度。不理这个差,扩散就漏,解就糊。

把 S_f 拆成两块#

面积矢量拆成两份。一份沿 $\mathbf{d}$,一份是余量。

$\mathbf{E}_f$ 是沿 $\mathbf{d}$ 的部分。$\mathbf{T}_f$ 是切向余量。代进通量,项就分成两截。

前一项只用 $\phi_N-\phi_P$。隐式进矩阵。只有两个单元中心。后一项要面平均梯度 $\overline{\nabla\phi}_f$,由相邻梯度平均得到。很难写进矩阵系数,于是作为显式修正扔到右端。用最小二乘把这个梯度重新拼起来的做法,写在在非结构网格上把梯度留住的移动最小二乘里。延迟修正发抖时,犯人有时不是 ρ,而是这个梯度。

$\theta=0$ 时 $\mathbf{T}_f=0$。退回原来的中心差分。

最小、正交、过松弛：三种分解

$\mathbf{E}_f$ 沿 $\mathbf{d}$ 的长度不是唯一的。记 $\hat{\mathbf{d}}=\mathbf{d}/|\mathbf{d}|$, $\cos\theta=\hat{\mathbf{d}}\cdot\mathbf{S}_f/|\mathbf{S}_f|$,有三种。

• 最小修正: $|\mathbf{E}_f|=|\mathbf{S}_f|\cos\theta$。$\mathbf{T}_f$ 最短。
• 正交修正: $|\mathbf{E}_f|=|\mathbf{S}_f|$。隐式大小保持为面积。
• 过松弛: $|\mathbf{E}_f|=|\mathbf{S}_f|/\cos\theta$。非正交越大,隐式部分越被拉长。

上面的画面里把 θ 推到60°,过松弛的 $\mathbf{E}_f$ 会长过面积矢量。对角系数变大了。最小修正相反,$\mathbf{E}_f$ 变短,$\mathbf{T}_f$ 变长。

代码：延迟修正拉普拉斯求解器

显式修正用的是上一次迭代的梯度。解一次不算完。修正要更新着扫很多遍。这就是延迟修正。

import numpy as np
 
def decompose_face(S_f, d, scheme="over"):
    """面积矢量 S_f = E_f（隐式，沿 d） + T_f（显式修正）。"""
    d_hat = d / np.linalg.norm(d)
    Smag = np.linalg.norm(S_f)
    dS = d_hat @ S_f                      # = |S_f| cos(非正交角)
    if scheme == "min":                   # 最小修正
        E = dS * d_hat
    elif scheme == "orth":                # 正交修正
        E = Smag * d_hat
    else:                                 # 过松弛 (over-relaxed)
        E = (Smag**2 / dS) * d_hat
    return E, S_f - E                     # E_f, T_f
 
def sweep_residual(S_f, d, scheme, sweeps=25):
    """E_f 隐式，T_f 延迟。残差以 rho = |T|/|E| 的比率收缩。"""
    E, T = decompose_face(S_f, d, scheme)
    rho = np.linalg.norm(T) / np.linalg.norm(E)
    r, hist = 1.0, [1.0]
    for _ in range(sweeps):
        r *= rho                          # 一次延迟修正 = 残差 x rho
        hist.append(r)
    return rho, hist
 
S_f = np.array([1.0, 0.0])                # 面法向 (面积矢量, |S|=1)
for ang in (20, 45, 65):
    d = np.array([np.cos(np.radians(ang)), np.sin(np.radians(ang))])
    for s in ("min", "orth", "over"):
        rho, hist = sweep_residual(S_f, d, s)
        tag = "收敛" if rho < 1 else "发散"
        print(f"theta={ang:2d}  {s:5s}  rho={rho:.3f}  {tag}  (25 次后 {hist[-1]:.1e})")

输出的意思是,角一大,三种分法就分开。20°时三个都还温和。65°时最小修正以 $\rho=\tan 65°\approx 2.14$ 发散,过松弛仍以 $\rho=\sin 65°\approx 0.91$ 收缩。同一张网格,同一个方程,拆法决定收敛还是发散。

三种 ρ 只由角度闭合。面积长度取 1 时,最小修正里 E 的长度是 cosθ,T 的长度是 sinθ,所以 ρ 是 tanθ。tan 45° 等于 1。45°就是收敛的墙。在这个角上,延迟循环给残差乘 1。扫 25 次,数字原样不动。20°的 tan 约为 0.364。65°的 tan 约为 2.14。发散一侧的数字来自这个正切。

正交修正里 E 的长度是 1,T 的长度是 2 减去 2cosθ 之后的平方根。ρ 要等于 1,就得 cosθ 为 1/2,也就是 60°。65°时已经约为 1.07。要等到 70°才说正交危险,说晚了。墙在 60°。

过松弛的 ρ 是 sinθ。65°时约为 0.906。代码片段把残差乘这个比值 25 次。0.906 的 25 次方约为 0.085。25 次之后显式误差大约还剩 8%。同一式子用在 20°,sin 约为 0.342,25 次方掉到 $10^{-12}$ 量级。把修正次数习惯地设成 2 或 3,在 65°的面上是 0.906 的立方,约为 0.74。大约四分之三的误差留在这个时间步里。

代码里的角是 d 与 S 之间的角。把 S 固定成单位矢量,20°、45°、65°就是比较点。最小修正在 45°的正切恰好是 1,不是凭经验挑的门槛。

为什么 ρ 是 T 的长度除以 E 的长度。延迟修正按 E 的长度做隐式系数,每次扫描把与 T 的长度成比例的误差显式送回去。每扫一次,残差按这个比值缩或胀。比值小于 1,循环才关得上。

非正交变大时会怎样

三种残差放在同一画面里一起走。角度滑块改 ρ,每扫一次曲线就按这个比值下降。

30°以下三条一起掉。45°时只有最小修正变成水平。过了 60°,正交的 ρ 也越过 1,只剩过松弛还停在 sinθ 上。

代价是有的。过松弛把隐式部分拉长,显式修正也跟着大。要碰到同样的残差,得多扫几遍。nNonOrthogonalCorrectors 就是这个次数。网格好,0 或 1 就够。歪得厉害,会被推到 2 或 3。65°上这两三次仍然留下大约四分之三的误差。一味加大,一个时间步就变贵。先把角本身压小,几乎总是更便宜。

为何在45°放弃最小修正#

45°上最小修正不再减小残差,因为 ρ 恰好是 1。正交修正在 60°撞上同一堵墙,65°时已经越过去了。过松弛的 ρ 是 sinθ,直到接近 90°都不会超过 1。一步之内留下的误差,就是这笔代价。

是再加扫描次数,还是在角超过 60°之前重铺网格。三角形覆盖点的顺序决定这个角。用两种方式把点云铺成三角形问的是,这个顺序在哪里留下非正交。加长修正循环更便宜,还是把前沿再推一遍更便宜。