网格歪斜时梯度也跟着偏 —— Green–Gauss 与加权最小二乘重构

线性场在乱模板上仍能对到小数点后六位。曲面则要靠n=1压住远处邻居，误差才会下来。

残差先是诡异地趴平，随后发散了。翻日志才发现，靠近壁面的一个单元里，温度梯度和真实方向几乎成了直角。再看网格，就那个单元格外扁。数值本身没问题。是构造梯度的方法被网格形状绊倒了。

在有限体积法（FVM）里，单元梯度是扩散通量、二阶重构和限制器的输入。仅靠单元中心值不够。要汇集邻居信息才能估计 $\nabla\phi_C$ 。拿来对照的是 Green–Gauss 和最小二乘。畸变网格上谁先歪、距离加权把误差拉回多少，都对着输出里的数字看。

把 anisotropy 降到 0.2 附近，模板会被压扁。把 weight exponent 从 0 提到 1，橙色箭头会转回绿色真值。

angle error 16.8° · |∇φ| 0.64 vs 1.08 true — the reconstructed gradient is skewed

Stencil anisotropy: 0.35Field curvature κ: 1.40Weight exponent n: 0.0Neighbor clustering: 0.60

同时把 curvature κ 抬高，两支箭头的夹角会张开。扁模板上，只看加权能不能把方向拧回来。

Green–Gauss：用包围的面来平均#

把散度定理直接用到单元上。将体积分换成面积分。

\nabla \phi_C = \frac{1}{V_C} \sum_{f} \phi_f \, \mathbf{n}_f \, A_f

$V_C$ 是单元体积， $\mathbf{n}_f$ 是面外法线， $A_f$ 是面积， $\phi_f$ 是面上的值。

面值 $\phi_f$ 无法直接得到。要从单元中心值插值。最常见的线性插值是这样：

\phi_f = g_f\, \phi_C + (1 - g_f)\, \phi_N

$g_f$ 是按面到邻居中心 $N$ 的距离比确定的权重。计算便宜且直观。在结构网格上能达到二阶精度。

畸变网格动摇了 Green–Gauss#

问题藏在 $g_f$ 里。线性插值假设面中心和两个单元中心在同一条直线上。在非结构网格里，这个假设几乎总是被打破。

如果面中心偏离连接两个单元中心的连线（ $\text{skewness}$ ，歪斜），插值得到的 $\phi_f$ 就有偏差。单元越扁、邻居尺寸差别越大，偏差越大。这个误差在迭代中累积。最终梯度指向错误的方向。

加上修正项能缓解。但修正又需要梯度。于是出现循环。迭代修正代价高，畸变严重时收敛变慢。

最小二乘：给邻居拟合一个平面

换个思路。不经过面。直接对单元 $C$ 周围的值分布拟合一个平面。

若邻居 $k$ 位于偏移量 $\mathbf{d}_k = \mathbf{x}_k - \mathbf{x}_C$ 处，线性近似为 $\phi_k \approx \phi_C + \nabla\phi_C \cdot \mathbf{d}_k$ 。对所有邻居，最小化这个误差的平方和。

\min_{\nabla\phi_C} \sum_k w_k \left( \phi_k - \phi_C - \nabla\phi_C \cdot \mathbf{d}_k \right)^2

$w_k$ 是各邻居的权重。求导令其为零，得到正规方程。

\left( \sum_k w_k\, \mathbf{d}_k \mathbf{d}_k^{\mathsf T} \right) \nabla\phi_C = \sum_k w_k \left( \phi_k - \phi_C \right) \mathbf{d}_k

二维时是一个 $2\times2$ 系统。展开写出来是这样：

\begin{pmatrix} \sum w\, d_x^2 & \sum w\, d_x d_y \\[2pt] \sum w\, d_x d_y & \sum w\, d_y^2 \end{pmatrix} \begin{pmatrix} \partial_x\phi \\[2pt] \partial_y\phi \end{pmatrix} = \begin{pmatrix} \sum w\, \Delta\phi\, d_x \\[2pt] \sum w\, \Delta\phi\, d_y \end{pmatrix}

$\Delta\phi = \phi_k - \phi_C$ 。左边矩阵只由邻居布局决定。可以每个单元算一次并存下来。优点很大。没有面插值假设。对畸变不敏感。只要邻居有 3 个以上就总能求解。

面插值垮掉的地方，和歪斜单元上的 Green–Gauss 与最小二乘是同一处。邻居只有三个也能解正规方程，可这三个点若排成一条线，横向的答案就是空的。

用距离加权

权重 $w_k$ 决定精度。常见的选择是距离的倒数。

w_k = \frac{1}{\lvert \mathbf{d}_k \rvert^{\,n}}

$n=0$ 时对所有邻居一视同仁。远处的邻居值差大，相应的曲率误差也大。那个邻居会把平面拟合往自己那边拉。把 $n$ 提到 1，远处邻居的话语权就减小。权重落在局部线性假设更成立的近邻上。无粘问题通常用 $n=1$ 。这一档压低远处邻居的表决权，曲率误差跟着下降。

Python —— 线性场精确，曲面靠加权#

分两步验证最小二乘的威力。先在随意撒开的模板上重构线性场 $\phi = 1.3x - 0.7y + 2$ 的梯度。理论上误差应为 0。接着在曲面 $\phi = e^{0.7x}\cos(0.9y)$ 上，让无加权一方与 $n=1$ 加权一较高下。

import numpy as np
 
def weighted_lsq_gradient(offsets, dphi, n):
    # offsets: (K,2) 邻居偏移量 d_k；dphi: (K,) 差值 phi_k - phi_C
    w = 1.0 / np.maximum(np.linalg.norm(offsets, axis=1), 1e-12) ** n
    M = (offsets[:, :, None] * offsets[:, None, :] * w[:, None, None]).sum(axis=0)
    b = (offsets * (dphi * w)[:, None]).sum(axis=0)
    return np.linalg.solve(M, b)
 
xc = np.array([0.4, -0.2])
 
# (A) 线性场 —— 无论何种模板，最小二乘都精确
lin = lambda x, y: 1.30 * x - 0.70 * y + 2.0
scatter = np.array([[0.45, 0.30], [-0.40, 0.12], [0.05, -0.48], [-0.33, -0.29],
                    [0.90, 0.15], [-0.11, 0.44], [0.28, -0.09]])   # 不规则模板
dlin = np.array([lin(*(xc + d)) - lin(*xc) for d in scatter])
gA = weighted_lsq_gradient(scatter, dlin, 1.0)
print(f"linear field   LSQ = ({gA[0]:+.6f}, {gA[1]:+.6f})   exact = (+1.300000, -0.700000)")
 
# (B) 曲面 —— 距离加权削减重构误差
phi  = lambda x, y: np.exp(0.7*x) * np.cos(0.9*y)
grad = lambda x, y: np.array([0.7*np.exp(0.7*x)*np.cos(0.9*y),
                             -0.9*np.exp(0.7*x)*np.sin(0.9*y)])
near = np.array([[0.10, 0.05], [-0.08, 0.09], [0.06, -0.10], [-0.09, -0.06], [0.11, 0.10]])
offs = np.vstack([near, [0.90, 0.70]])          # 一个远处的邻居
dphi = np.array([phi(*(xc + d)) - phi(*xc) for d in offs])
gt   = grad(*xc)
for n in (0.0, 1.0):
    g = weighted_lsq_gradient(offs, dphi, n)
    err = np.linalg.norm(g - gt)
    print(f"curved field   n={n:.0f}  LSQ = ({g[0]:+.3f}, {g[1]:+.3f})   |error| = {err:.4f}")
print(f"curved field   true = ({gt[0]:+.3f}, {gt[1]:+.3f})")

输出是这样的。

linear field   LSQ = (+1.300000, -0.700000)   exact = (+1.300000, -0.700000)
curved field   n=0  LSQ = (+0.885, +0.200)   |error| = 0.0293
curved field   n=1  LSQ = (+0.891, +0.205)   |error| = 0.0222
curved field   true = (+0.911, +0.213)

线性场上，模板再扭，梯度也精确到小数点后六位。输出里的 (+1.300000, -0.700000) 就是这个值。线性函数在每个邻居上都让差值和梯度的点积严格成立，所以不规则模板里的 (0.90, 0.15) 并不制造误差。畸变在这个场上推不动解。

曲面就不一样。真值是 (+0.911, +0.213)。n=0 的 (+0.885, +0.200) 在 x 上少 0.026，在 y 上少 0.013。n=1 的 (+0.891, +0.205) 只少 0.020 和 0.008。范数从 0.0293 降到 0.0222。差值 0.0071，大约 24%，和四分之一是同一件事。远处邻居 (0.90, 0.70) 的分量比近处五点的 0.1 量级大得多。n=0 时它和别人投同样的票，把平面拽走。n=1 时距离把这票削掉。

邻居排成一条线时，别信那张逆矩阵

最小二乘抗畸变，但不是万能的。左边矩阵 $M = \sum w_k \mathbf{d}_k \mathbf{d}_k^{\mathsf T}$ 可能病态。邻居排成一条直线时， $M$ 接近奇异。沿模板方向梯度定得好，横向几乎是自由的。条件数会爆。

把 Collinearity 提到 0.9 附近，邻居就挤成一列，cond(M) 从几十冲到无穷大。

cond(M) = 1.4 (λ_max 3.14 / λ_min 2.271) — well-conditioned

Collinearity: 0.20Neighbor count: 6

椭圆被拉得又长又细的那根轴，就是梯度定不住的轴。贴壁的边界层单元正是这个形状。对角元是偏移平方之和，邻居只沿一个方向拉长，另一条对角就空了。边界单元的邻居堆在流体一侧。把壁对面的邻居或面心补进模板，空轴会稍微填上，条件数下来。网格固定时， $M$ 不随速度场改变。不必每步重装，把逆矩阵或 QR 存下来。距离指数 $n$ 看问题。粘性项大的时候， $n=2$ 有时更合适。

移动最小二乘是把这份权重按邻居重做的下一步。Green–Gauss 便宜，但怕偏斜。加权最小二乘不怕偏斜，可模板一旦排成一条线，矩阵会先病掉。选方法之前，先看这个单元的邻居是不是铺在一个平面上。