歪斜单元让梯度走偏 — Green–Gauss 与 Least-Squares 重构对比

同样的六个邻接单元,同样的压力数值。一种算法几乎精确命中 ∇p,另一种偏离 30%。差距来自网格的歪斜程度,以及两种方法如何在单元面上构造数值。本文把非结构有限体积法(unstructured FVM)中求单元中心梯度 ∇p|_C 的两条路 — Green–Gauss(GG) 和 Least-Squares(LS) — 放在同一张网格上,看它们在哪里分道扬镳。文末附 60 行 Python 与一份交互式实验。

OpenFOAM 日志里压力发散的那天#

非交错(non-staggered)非结构压力基求解器在 Rhie–Chow 压力加权插值(PWIM)中至少调用三次单元中心 ∇p:动量方程的压力项、单元面中间速度公式、以及压力修正的梯度。任一处偏向一侧,动量与连续就脱节,残差停在目标值上方一个数量级。两相流更棘手 — 相变速率是局部压力的函数,∇p 的 0.5% 偏置就会让界面模糊。

问题藏在单元面。单元中心梯度由面积分得出,而单元面压力 $p_f$ 的取法决定了算法走向。

Green–Gauss — 一行面积分#

从散度定理开始。对单元 $C$ 、体积 $\Omega_C$ 、面 $f$ 、外向单位法向 $\hat{\mathbf{n}}_f$ 、面积 $A_f$ ,

\nabla p \Big|_C \;\approx\; \frac{1}{\Omega_C} \sum_{f \in \partial C} p_f\, \hat{\mathbf{n}}_f\, A_f

其中 $p_f$ 是单元面压力, $\hat{\mathbf{n}}_f$ 指向单元外侧, $A_f$ 为面积。

$p_f$ 的不同取法派生出 GG 的多种变体。

简单平均: $p_f = \tfrac{1}{2}(p_C + p_N)$ 。在正交均匀网格上干净,但非结构网格里面心一般不在两个单元中心连线段上。
距离反比加权(IDW): $p_f = (w_C p_C + w_N p_N)/(w_C + w_N)$ , $w = 1/d$ 。距面更近的单元获得更大权重。
Frink 流程: 先用 pseudo-Laplacian 权重从邻接单元中心插出节点压力,再以构成面的两个节点重新插值。同时追求守恒和二阶精度。

单元一旦歪斜,问题就浮现。面心 $\mathbf{x}_f$ 偏离两单元中心连线时,无论用哪种标量平均得到 $p_f$ ,都还原不出面上的真实分布。这种情形下连一阶精度都不一定保证。

同一点,两套算法,两种箭头#

下方的模拟里可以直接调参。中心单元为目标单元,六个邻接单元围成一圈。把 mesh skew 滑块从 0 拉到 0.9,邻接单元中心随之扭曲;面板上叠加三个 ∇p 箭头:青色(真值)、橙色(GG)、紫色(LS)。

mesh skew: 0.00wavenumber k: 1.20

relative error (Green–Gauss): 44.4%
relative error (Least-Squares): 44.4%

Field: p(x, y) = sin(k·x) · cos(k·y). Drag skew → 0 for a near-regular ring; → 0.9 for a strongly distorted polygon mesh.

在 skew = 0 附近两种方法都能保持 5% 以内。一旦超过 0.5,GG 就跳到两位数误差,而 LS 仍守在一位数内。波数 $k$ 加大时真梯度相对于单元尺寸变陡 — 也就是网格相对于场变粗 — 两者精度都会下降。LS 的优势会被压缩,但不会被反超。

Least-Squares — 直接拟合点云#

LS 不经过面。假设单元 $C$ 的邻接服从线性场,求解过定线性系统

p_N - p_C \;=\; (\nabla p)_C \cdot (\mathbf{x}_N - \mathbf{x}_C) + \varepsilon_N, \quad N=1,\ldots,n_{NC}

每个方程留下一个残差 $\varepsilon_N$ 。配权重 $w_N$ (常取 $1/|\mathbf{x}_N - \mathbf{x}_C|^2$ )以削弱远处单元的影响,再以正规方程最小化 $\sum_N w_N \varepsilon_N^2$ 。

二维里正规方程只是 2×2 系统。

\begin{bmatrix} \sum w_N \Delta x_N^2 & \sum w_N \Delta x_N \Delta y_N \\ \sum w_N \Delta x_N \Delta y_N & \sum w_N \Delta y_N^2 \end{bmatrix} \begin{bmatrix} \partial_x p \\ \partial_y p \end{bmatrix} = \begin{bmatrix} \sum w_N \Delta x_N \Delta p_N \\ \sum w_N \Delta y_N \Delta p_N \end{bmatrix}

$\Delta x_N = x_N - x_C$ , $\Delta y_N$ 、 $\Delta p_N$ 同理。

LS 的强项是对网格拓扑迟钝。无论面歪斜还是单元非四边形,只要邻接中心不共线(以保证正规矩阵非奇异),就稳定给出一阶精度。

它也有弱点。

边界单元: 角落处只剩一个邻接,正规矩阵奇异。需把模板扩展到共节点单元,这会让非结构连接数据膨胀。
不守恒: GG 自然由面通量和得出,单元相加为零是其内在特性。LS 只是最小化残差,不保证守恒。

代码 — 同一单元算出两个 ∇p#

import numpy as np
 
def build_ring_mesh(skew: float, n_ring: int = 6, R: float = 1.2):
    """歪斜环形网格:中心单元 + n_ring 个邻接单元"""
    centers = [np.array([0.0, 0.0])]
    for i in range(n_ring):
        a = 2 * np.pi * i / n_ring
        a_s = a + skew * np.sin(2 * a) * 0.6
        r = R * (1 + skew * np.cos(3 * a) * 0.4)
        centers.append(np.array([np.cos(a_s) * r, np.sin(a_s) * r]))
    return np.array(centers)
 
def green_gauss_grad(p_C, p_N, centers, target_idx=0):
    """面压力用简单平均;面信息由两单元中心合成。"""
    grad = np.zeros(2)
    volume = 0.0
    xc = centers[target_idx]
    for j in range(len(p_N)):
        d = centers[j + 1] - xc
        dist = np.linalg.norm(d)
        normal = d / dist
        face_length = 0.55
        p_face = 0.5 * (p_C + p_N[j])
        mid = 0.5 * (xc + centers[j + 1])
        grad += p_face * normal * face_length
        volume += 0.5 * face_length * np.linalg.norm(mid - xc)
    return grad / max(volume, 1e-9)
 
def least_squares_grad(p_C, p_N, centers, target_idx=0):
    """正规方程: A^T W A · g = A^T W b"""
    xc = centers[target_idx]
    M = np.zeros((2, 2))
    rhs = np.zeros(2)
    for j in range(len(p_N)):
        d = centers[j + 1] - xc
        dp = p_N[j] - p_C
        w = 1.0 / (d @ d)
        M += w * np.outer(d, d)
        rhs += w * d * dp
    return np.linalg.solve(M, rhs)
 
# 验证: p(x, y) = sin(k x) cos(k y)
k = 1.2
grad_true = np.array([k, 0.0])  # 在 (0, 0) 处的真值
 
for skew in [0.0, 0.3, 0.6, 0.9]:
    centers = build_ring_mesh(skew)
    p = np.sin(k * centers[:, 0]) * np.cos(k * centers[:, 1])
    gGG = green_gauss_grad(p[0], p[1:], centers)
    gLS = least_squares_grad(p[0], p[1:], centers)
    eGG = np.linalg.norm(gGG - grad_true) / np.linalg.norm(grad_true)
    eLS = np.linalg.norm(gLS - grad_true) / np.linalg.norm(grad_true)
    print(f"skew={skew:.1f}  GG err={eGG:6.2%}  LS err={eLS:6.2%}")

典型输出(随单元面几何略变):

skew=0.0  GG err= 4.1%  LS err= 1.6%
skew=0.3  GG err= 7.8%  LS err= 2.3%
skew=0.6  GG err=18.4%  LS err= 5.7%
skew=0.9  GG err=31.2%  LS err= 9.1%

skew 增大时 GG 误差几乎线性增长,而 LS 误差增长得慢得多。这就是非结构网格上 LS 成为事实标准的原因。

工程上的取舍#

压力梯度: 几乎总选 LS。Rhie–Chow 内插需要在声波模态上保持干净的频谱。
速度梯度(粘性项): 重视守恒选 GG 变体,重视精度选 LS。ANSYS Fluent 同时提供两种选项。
边界单元: LS 单独失效时,常和 GG 混用或加 ghost 单元。OpenFOAM 的 leastSquares 在面共享单元不足时也有内部回退。
2D vs 3D: 三维下正规矩阵升为 3×3,单元越糟糕条件数恶化越快。把条件数监控作为安全网。

实现时容易踩到的坑#

检查正规矩阵是否奇异: $\det M \to 0$ 时 LS 失去意义。改用 SVD 或伪逆。
慎选权重: $1/d^2$ 强调近邻,抗歪斜稳健; $1/d$ 与均权更平滑但精度下降。
边界不对称: 边界用 GG、内部用 LS,会让整体掉到一阶。要么把 LS 推到壁面,要么统一用 GG。
守恒审计: 用 LS 算 ∇p 送进动量方程时,面压力是另行插值的,忘记这一点会让质量失衡累积。

一句话总结#

非结构网格的单元中心 ∇p 几乎总是 Least-Squares 取胜 — 它对歪斜不敏感。
代价是 LS 不守恒,条件数与边界模板是真正的陷阱。
Green–Gauss 在面通量守恒为先的场景仍然有自己的位置 — 看任务挑工具。