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

在边界模板上切断最小二乘的那一刻，内部守住的个位数误差也连同阶数一起掉下去。

同样六个邻居，同样的压力。一支箭几乎指住 $\nabla p$ ，另一支偏出去不少。差别在于单元歪了多少，以及面上的压力是怎么造出来的。

把 mesh skew 滑块从0推到0.9，橙色的 Green–Gauss 箭头离开真值，紫色的最小二乘留得更短。

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时两者都在个位数。过了0.5，只有橙色升到两位数。表上的数字在代码后面。

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

非交错的非结构压力基求解器，在 Rhie–Chow 压力加权插值里至少三次调用单元中心的 $\nabla p$ 。动量里的压力项，面上的中间速度，压力修正的梯度。只要有一个歪向一侧，动量和连续就分开，残差停在高一个量级的地方。两相里相变率是局部压力的函数， $\nabla 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$ 的做法分开 Green–Gauss 的变体。

简单平均: $p_f=\tfrac12(p_C+p_N)$ 。正交均匀网格上干净。非结构里，面心并不在两个单元中心连线的中点。
距离反权: $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$ 用哪种平均，都抓不住面上的分布。一阶常常都保证不了。分母里的 $\Omega_C$ 若是假面积，守恒就只剩形式。用对角线外积把单元面积闭上的地方，就是这块面积要和 $\partial x/\partial x=1$ 咬合的地方。

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

前面格子的环是一个中心单元加六个邻居。skew 为0时两支箭贴着青色真值。过了0.5，橙色先离开。把波数 $k$ 加大，真 $\nabla p$ 相对单元尺寸变陡，两种方法都受伤。那时紫色的优势变窄。

Least-Squares — 直接拟合点云#

最小二乘不经过面。假定邻居满足线性分布，解超定方程组。

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 w_N\varepsilon_N^2$ 最小的 $\nabla p$ 来自正规方程。

二维正规方程是 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$ 是同一种差。

邻居中心不只落在一条直线上，正规矩阵就非奇异，面歪了也还能保住一阶。这就是它对拓扑迟钝的原因。

弱点也在。只有一个邻居的边界角点上，正规矩阵奇异。模板得扩到共享节点的单元，非结构上这个索引很重。而且没有守恒。Green–Gauss 是面通量之和，单元合起来容易得到0。最小二乘只是最小化残差，不保证那个和。

代码 — 同一单元算出两个 ∇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 为0时最小二乘的1.6%不是歪斜。代码里 $k=1.2$ ，环半径 $R=1.2$ ，于是 $kR=1.44$ rad。场是 $\sin(kx)\cos(ky)$ ，原点的真梯度是 $(k,0)=(1.2,0)$ 。线性拟合装不下这个曲率。skew 为0时正规矩阵只剩下角度。权重 $w=1/|\mathbf{d}|^2$ 时， $w\Delta x^2$ 的和就是 $\cos^2$ 的和。间隔60度的六个，这个和是3，交叉项是0，所以 $M=3I$ 。条件数是1。1.6%不是矩阵病了的结果。

表里的比是2.6、3.4、3.2、3.4。即 $4.1/1.6$ 、 $7.8/2.3$ 、 $18.4/5.7$ 、 $31.2/9.1$ 。skew 从0到0.9，Green–Gauss 变成 $31.2/4.1\approx 7.6$ 倍，最小二乘是 $9.1/1.6\approx 5.7$ 倍。Green–Gauss 的增量是 $+3.7$ 、 $+10.6$ 、 $+12.8$ ，不是直线。后面更陡。

skew 为0.9时，代码里的半径是 $R(1\pm 0.36)$ ，角度摆动的幅度是 $0.54$ rad。就是 $0.9\times 0.4$ 和 $0.9\times 0.6$ 。面长 $0.55$ 在代码里是固定的。六个面共用这个长度时，它在分子和分母里约掉。把0.55换成另一个常数，印出的4.1%和31%不动。这两个百分数属于环的角度和正弦场。

skew为0时，到邻居的距离是 $R=1.2$ ，面中点在一半处，0.6。去掉长度之后，每个面的分母贡献是 $0.5\times 0.6$ ，六面合计1.8。代码累加的体积0.99是再乘一次0.55，进不了商。邻居间距为 $R$ 的正六边形面积是 $(\sqrt{3}/2)R^2$ 。 $R^2=1.44$ 时大约1.25。约掉之后的分母1.8比这块面积大44%。若分子也等于真正的面积分，梯度会落到真值的 $1/1.44$ ，相对误差是 $1-1/1.44\approx 0.31$ 。印出的 Green–Gauss 误差是4.1%，比这个0.31小得多。分母偏大的同一比例上，面压力之和也偏大，两份几何误差抵掉大部分之后，才剩下4.1%。最小二乘的1.6%与这块面积无关。

弯曲的大小可以从 $k$ 和 $R$ 读出来。 $kR=1.44$ ，泰勒三次项 $(kR)^3/6$ 大约0.50，约是线性项 $kR$ 的三分之一。环上的样本已经弯了这么多。线性拟合留下的梯度误差只有1.6%，这点弯曲的大部分在对称模板上离开了梯度。skew 一升高，这份抵消就破了。从4.1%到31.2%的7.6倍不是新长出来的曲率。那是本来抵掉的几何误差，随着角度和半径的摆动露了出来。

工程上的取舍

压力梯度几乎总是最小二乘。Rhie–Chow 里要把声波模式抓干净。

速度梯度，也就是粘性项，是分开的。要守恒就用 Green–Gauss 的变体，精度优先就用最小二乘。Fluent 两个都放。

边界单元若单独解不了最小二乘，就用 Green–Gauss 补，或者放 ghost。OpenFOAM 的 leastSquares 在共享面的单元不够时，内部也有 fallback。

三维里正规矩阵是 3×3，单元越糙，条件数坏得越快。看着这个数，就是保险。

实现时容易踩到的坑

$\det M$ 贴到0，最小二乘就没有意义。改成 SVD 或伪逆。

权重 $1/d^2$ 放大近处单元，对歪斜迟钝。 $1/d$ 或均匀权重更光滑，精度更低。

只把边界落到 Green–Gauss、内部留最小二乘，整个系统会掉到一阶。要么把最小二乘推到边界，要么全部统一成 Green–Gauss。

即使用最小二乘解了压力梯度，面压力仍然另做插值。忘掉这一点，质量不平衡会堆起来。

在边界切断最小二乘，9%也一起塌#

skew 为0.9的那一行把选择变成数字。选守恒，就吃下 Green–Gauss 的31%。选精度，就留下最小二乘的9%。可是在边界切断这个最小二乘的一刻，内部守住的个位数会连同阶数一起沉下去。条件数为1的正六边形环，在边界角点上并不成立。

梯度箭头对了，面速度若记住了松弛系数，收敛解本身也会变。日志发散时，剩下的问题是箭头错了，还是面记住了松弛。Rhie–Chow 记住欠松弛的地方就是这个问题。