把矩形网格裹到机翼上 — 曲线坐标变换与网格度量

把 Fhat 直接写成 y_η 和 x_η，混合导数会相消，残差落到舍入。先除以 J 再乘回去，这个相消就会错开。

一位工程师想求解机翼周围的流动。手上的求解器只认识直角网格。可机翼是曲面。无论把矩形单元切得多细，曲面也只能画成阶梯状。出路不是把网格弯过去，而是把空间本身弯过去。把扭曲的物理网格映射回一张笔直的逻辑网格，再在那张笔直的网格上求解方程。本文跟随让这次往返成为可能的坐标变换与网格度量。读完之后，你会明白雅可比为何能捕捉网格折叠，以及曲线网格上一股均匀流要保持均匀需要什么——并用numpy直接验证。

在网格映射器里提高 shear 和 amplitude，单元会扭起来。Jacobian heatmap 按 $J/J_0$ 给每一格上色。

shear (η 기울임) 0.30

amplitude 0.15

wave number 2

grid n = 16

J/J₀ min = 1.00, max = 1.00 · 접힘 없음 (J>0)

amplitude 收到 0.3 附近、wave number 收到 4 以上，红色单元出现。那些格子里的 $J$ 是负的。生成器越过这条线，度量的分母就会穿过 0。

物理域与逻辑域

曲线坐标法从两个域出发。物理域是包裹机翼的扭曲网格所在之处，坐标记作 $(x, y)$ 。逻辑域是把那张网格摊平后的单位正方形，坐标记作 $(\xi, \eta)$ 。

映射是连接两者的函数。

x = x(\xi, \eta), \qquad y = y(\xi, \eta)

其中 $\xi$ 对应网格某一方向的索引， $\eta$ 对应另一方向。在逻辑域中每个单元都是边长 $\Delta\xi = \Delta\eta = 1/N$ 的正方形。求解器只在这张笔直的网格上计算差分。所有扭曲信息都被吸收进映射函数里。

度量:网格建立的汇率表

麻烦在于，我们要解的方程是用关于 $x, y$ 的导数写成的。网格却沿 $\xi, \eta$ 移动。连接两者的汇率表就是网格度量(metric)。

从链式法则开始。

\frac{\partial}{\partial x} = \xi_x \frac{\partial}{\partial \xi} + \eta_x \frac{\partial}{\partial \eta}, \qquad \frac{\partial}{\partial y} = \xi_y \frac{\partial}{\partial \xi} + \eta_y \frac{\partial}{\partial \eta}

像 $\xi_x = \partial\xi/\partial x$ 这样的项就是度量系数。但我们手上握着的是反方向，即 $x_\xi, x_\eta, y_\xi, y_\eta$ ，因为把映射对 $\xi, \eta$ 求导即可直接得到。两者通过逆矩阵关系绑定。

\xi_x = \frac{y_\eta}{J}, \quad \xi_y = -\frac{x_\eta}{J}, \quad \eta_x = -\frac{y_\xi}{J}, \quad \eta_y = \frac{x_\xi}{J}

每个符号都是映射的一阶导数， $J$ 是下一节的雅可比。归纳起来，度量就是"容易得到的 $x_\xi$ 一族除以雅可比"。

雅可比与网格折叠

雅可比行列式衡量映射把面积放大或缩小了多少。

J = \det \begin{pmatrix} x_\xi & x_\eta \\ y_\xi & y_\eta \end{pmatrix} = x_\xi y_\eta - x_\eta y_\xi

其中 $J$ 是逻辑域中单位面积在物理域所占面积之比。 $J > 0$ 时映射保持定向。 $J = 0$ 时单元塌缩成一点或一线。 $J < 0$ 时单元被翻转。这就是网格折叠(grid folding)。在折叠网格上度量的分母穿过零而发散，求解器吐出 NaN。

$J = x_\xi y_\eta - x_\eta y_\xi$ 是单元两条边的叉积。和用对角线叉积量单元面积是同一项。上面 heatmap 变红的格子，就是 $J < 0$ 。

变换后的守恒形方程

现在把方程搬进逻辑域。看一个二维守恒形方程。

\frac{\partial Q}{\partial t} + \frac{\partial F}{\partial x} + \frac{\partial G}{\partial y} = 0

$Q$ 是守恒量， $F, G$ 是 $x, y$ 方向的通量。代入度量并乘以 $J$ 整理后得到强守恒形(strong conservation form)。

\frac{\partial (JQ)}{\partial t} + \frac{\partial \hat F}{\partial \xi} + \frac{\partial \hat G}{\partial \eta} = 0

变换后的通量为

\hat F = J(\xi_x F + \xi_y G), \qquad \hat G = J(\eta_x F + \eta_y G)

关键在于，像 $J\xi_x = y_\eta$ 、 $J\xi_y = -x_\eta$ 这样度量与雅可比之积会折回映射的简单导数。因此变换后的通量无需先除以度量再乘回。守住这一形式，曲线网格上质量与动量也能精确守恒。

用代码计算并验证度量

把理论搬到 numpy。定义一张波动的物理网格，用中心差分还原度量与雅可比，再看均匀通量是否真的守恒。

import numpy as np
 
def body_fitted_map(xi, eta, amp=0.15, shear=0.3, wavek=2.0):
    """把逻辑坐标 (xi, eta) ∈ [0,1]^2 映射到物理坐标 (x, y)。"""
    x = xi + shear * eta + amp * np.sin(np.pi * wavek * eta)
    y = eta + amp * np.sin(np.pi * wavek * xi)
    return x, y
 
def compute_metrics(x, y, dxi, deta):
    """用中心差分还原雅可比 J 与逆度量。"""
    x_xi  = np.gradient(x, dxi,  axis=1)   # ∂x/∂ξ
    x_eta = np.gradient(x, deta, axis=0)   # ∂x/∂η
    y_xi  = np.gradient(y, dxi,  axis=1)
    y_eta = np.gradient(y, deta, axis=0)
    J = x_xi * y_eta - x_eta * y_xi        # 雅可比行列式
    xi_x, xi_y  =  y_eta / J, -x_eta / J
    eta_x, eta_y = -y_xi / J,  x_xi / J
    return J, (xi_x, xi_y, eta_x, eta_y)
 
# 逻辑网格
N = 41
xi  = np.linspace(0, 1, N)
eta = np.linspace(0, 1, N)
dxi, deta = xi[1] - xi[0], eta[1] - eta[0]
XI, ETA = np.meshgrid(xi, eta)             # axis0=eta, axis1=xi
 
X, Y = body_fitted_map(XI, ETA, amp=0.15, shear=0.3, wavek=2.0)
J, _ = compute_metrics(X, Y, dxi, deta)
 
print(f"J 最小/最大 = {J.min():.4f} / {J.max():.4f}")
print("网格:", "发生折叠 (J<=0)" if J.min() <= 0 else "正常 (J>0)")

接着是自由来流守恒(freestream preservation)检验。代入常通量 $F, G$ ，残差应为零。把变换通量直接写成 $J\xi_x = y_\eta$ 的形式，中心差分的混合导数互相交换，残差落到舍入误差量级。

def freestream_residual(x, y, dxi, deta, F=1.0, G=0.5):
    """常通量 (F,G) 的自由来流残差，应为零。"""
    x_xi  = np.gradient(x, dxi,  axis=1)
    x_eta = np.gradient(x, deta, axis=0)
    y_xi  = np.gradient(y, dxi,  axis=1)
    y_eta = np.gradient(y, deta, axis=0)
    Fhat =  F * y_eta - G * x_eta          # = J(ξ_x F + ξ_y G)
    Ghat = -F * y_xi  + G * x_xi           # = J(η_x F + η_y G)
    dFhat = np.gradient(Fhat, dxi,  axis=1)
    dGhat = np.gradient(Ghat, deta, axis=0)
    return dFhat + dGhat
 
R = freestream_residual(X, Y, dxi, deta)
print(f"自由来流残差 max|R| (内部) = {np.abs(R[1:-1, 1:-1]).max():.2e}")
# 例) J 最小/最大 = 0.62 / 1.38, 自由来流残差 max|R| ≈ 1e-16

常见错误是用非守恒形求解:先把度量除以 $J$ ，再乘回通量。三维中这种不匹配会变大，在均匀流里播下虚假源项。Thomas 与 Lombard(1979)提出的对称守恒形度量可以防止它。二维中，像上面那样把 $\hat F$ 直接写成 $y_\eta, x_\eta$ 就足够了。

先看代码注释里的数。amp=0.15、shear=0.3、wavek=2、N=41 时， $J$ 落在 0.62 到 1.38。最大是最小的 $1.38/0.62$ ，约 2.23 倍。单元面积起伏超过两倍， $J$ 仍然全是正的。没有折叠的格子。同一张网上自由来流残差 $\max|R|$ 约 1e-16。这是 $F=1$ 、 $G=0.5$ 时的值。

面积在晃，残差却停在舍入，顺着公式看原因。 $\hat F = F y_\eta - G x_\eta$ ， $\hat G = -F y_\xi + G x_\xi$ 。把 $\xi$ 方向差分和 $\eta$ 方向差分相加，只要中心差分交换混合导数， $y_\eta$ 对 $\xi$ 的导数就和 $y_\xi$ 对 $\eta$ 的导数抵消。 $F$ 和 $G$ 是常数，可以提出去。剩下的是映射两个混合导数之差。同一个算子，这个差就是 0。把 $\xi_x$ 除以 $J$ 再乘回 $J$ ，解析上回到 $y_\eta$ 。离散时除法和乘法可能碰到不同位置的 $J$ 。这道缝就是均匀流里的假源。二维里把 $\hat F$ 直接写成 $y_\eta$ 、 $x_\eta$ ，缝就不会打开。注释里的 1e-16 是相消留下来的位置。

$\beta$ 也从同一式子里拿出数字。壁面 $\eta=0$ 处 $dy/d\eta = \beta/(e^\beta-1)$ ， $\eta=1$ 处是 $\beta e^\beta/(e^\beta-1)$ 。两个斜率之比是 $e^\beta$ 。 $\beta=4$ 时约 54.6 倍。正文说的数十倍就是这个比。 $\beta$ 趋于 0 时斜率趋于 1，也就是均匀网格。

在壁面附近聚集单元

要解边界层，就得把单元密集地聚在壁面附近。这也是同一坐标变换的一维版本。用指数函数拉伸均匀的 $\eta$ 来构造物理坐标 $y$ 。

y(\eta) = \frac{e^{\beta \eta} - 1}{e^{\beta} - 1}

其中 $\beta$ 是聚集系数。一维雅可比为 $dy/d\eta = \beta e^{\beta\eta}/(e^\beta - 1)$ ，该值小的地方网格变密。

在下面的模拟中调节 $\beta$ 。

clustering β = 2.5

nodes n = 16

첫 간격 Δy₀ = 0.0000 · 마지막 간격 = 0.0000 · 비율 = —×

把 $\beta$ 从 0 提到 4，壁面( $y=0$ )附近的第一个间距急剧收缩，与最后一个间距之比拉开到数十倍。这意味着同样的单元数能把边界层捕捉得好得多。

先除以J再乘回去，均匀流里会长出源项#

$J$ 为负的格子不接收通量。默认波浪，amp=0.15、shear=0.3、wavek=2、N=41， $J$ 在 0.62 和 1.38 之间。没有折叠。生成阶段只要出现一个红格，这张网就不进求解器。NaN 从那个分母里出来。

留下的形式只有一种。 $\hat F$ 用 $y_\eta$ 和 $x_\eta$ 来写。先除以 $J$ 再乘回去，二维也会裂开一道缝，三维里就是 Thomas–Lombard 挡住的假源。用同一套度量去转刚度本身，是另一件事，写在构成张量的坐标变换里。