第 5 章：网格和相应变换

本章解决的问题

CFD 最终需要在离散网格上求解控制方程。

对于简单矩形区域，可以直接采用 Cartesian 网格：

(xi,yj)(x_i,y_j)

但真实流体问题中的物体通常具有复杂曲面，例如：

  • 圆柱；
  • 翼型；
  • 叶片；
  • 船体。

如果仍然直接采用规则直角坐标网格，则物体表面可能不断切穿网格，使边界条件处理变得困难。

因此引入计算坐标，把复杂的物理区域，映射到规则的计算区域

之后在规则的计算空间网格上进行差分计算，再通过坐标变换解释其在物理空间中的含义。

网格变换的一般思想

简单起见以二维平面为例：

有正变换

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

以及反变换

ξ=ξ(x,y),η=η(x,y).\xi=\xi(x,y), \qquad \eta=\eta(x,y).

一个计算网格点 (ξi,ηj)(\xi_i,\eta_j)，对应物理空间中的一点(xij,yij)(x_{ij},y_{ij})

网格生成实际上是：

求出计算平面上所有规则网格点在物理平面中的坐标。

对于结构网格，即使物理空间中的网格线已经严重弯曲，只要仍然能够按照 (i,j)(i,j) 编号，并明确找到(i±1,j),(i,j±1)(i\pm1,j),\qquad(i,j\pm1) 这样的邻居，它仍然属于结构网格。

所以“结构网格”强调的是规则的拓扑连接关系。

微分变换

坐标变化以后，原本关于 x,yx,y 的导数必须重新写成关于 ξ,η\xi,\eta 的导数。

若

ϕ=ϕ(x,y)=Φ(ξ,η),\phi=\phi(x,y)=\Phi(\xi,\eta),

根据链式法则：

∂ϕ∂x=∂Φ∂ξ∂ξ∂x+∂Φ∂η∂η∂x,\frac{\partial\phi}{\partial x} = \frac{\partial\Phi}{\partial\xi} \frac{\partial\xi}{\partial x} + \frac{\partial\Phi}{\partial\eta} \frac{\partial\eta}{\partial x},

以及

∂ϕ∂y=∂Φ∂ξ∂ξ∂y+∂Φ∂η∂η∂y.\frac{\partial\phi}{\partial y} = \frac{\partial\Phi}{\partial\xi} \frac{\partial\xi}{\partial y} + \frac{\partial\Phi}{\partial\eta} \frac{\partial\eta}{\partial y}.

所以可以写成算子形式：

∂∂x=ξx∂∂ξ+ηx∂∂η\boxed{ \frac{\partial}{\partial x} = \xi_x\frac{\partial}{\partial\xi} + \eta_x\frac{\partial}{\partial\eta} } ∂∂y=ξy∂∂ξ+ηy∂∂η\boxed{ \frac{\partial}{\partial y} = \xi_y\frac{\partial}{\partial\xi} + \eta_y\frac{\partial}{\partial\eta} }

其中

ξx,ξy,ηx,ηy\xi_x,\xi_y,\eta_x,\eta_y

被称为坐标变换的度量系数。

变换矩阵和 Jacobian

定义从计算空间到物理空间的局部变换矩阵：

Js=[xξxηyξyη].\mathbf J_s= \begin{bmatrix} x_\xi & x_\eta\\ y_\xi & y_\eta \end{bmatrix}.

其行列式：

J=xξyη−xηyξ\boxed{ J=x_\xi y_\eta-x_\eta y_\xi }

称为 Jacobian。

这里采用的是

J=∂(x,y)∂(ξ,η)J= \frac{\partial(x,y)} {\partial(\xi,\eta)}

这一约定，不同教材也可能定义其倒数。

几何意义

Jacobian 的几何意义是：

计算平面中一个很小的面积单元，在映射到物理平面后被放大或缩小的倍数。

因此必须要求：

J≠0.J\neq0.

如果

J=0,J=0,

说明某个网格单元被压扁，坐标变换失去一一对应关系，网格出现退化。

逆变换系数

根据矩阵求逆：

[ξxξyηxηy]=[xξxηyξyη]−1.\begin{bmatrix} \xi_x & \xi_y\\ \eta_x & \eta_y \end{bmatrix} = \begin{bmatrix} x_\xi & x_\eta\\ y_\xi & y_\eta \end{bmatrix}^{-1}.

因此：

ξx=yηJ,ξy=−xηJ,\xi_x=\frac{y_\eta}{J}, \qquad \xi_y=-\frac{x_\eta}{J}, ηx=−yξJ,ηy=xξJ.\eta_x=-\frac{y_\xi}{J}, \qquad \eta_y=\frac{x_\xi}{J}.

实际计算中，通常更容易直接得到

x(ξ,η),y(ξ,η)x(\xi,\eta),\qquad y(\xi,\eta)

以及

xξ,xη,yξ,yη,x_\xi,x_\eta,y_\xi,y_\eta,

然后通过 Jacobian 得到逆变换系数。

极坐标贴体坐标示例

对于圆柱，可以选择：

ξ=θ,η=r.\xi=\theta,\qquad\eta=r.

于是：

x=ηcos⁡ξ,y=ηsin⁡ξ.x=\eta\cos\xi, \qquad y=\eta\sin\xi.

圆柱表面

r=Rr=R

直接变成计算空间中的：

η=R.\eta=R.

所以圆柱表面本身就是一条坐标线。

这说明：

极坐标本质上就是一种具有解析表达式的贴体网格。

一般翼型、船体等没有简单解析坐标系，所以才需要数值方法生成更加一般的曲线坐标。

强守恒形式

如果物理空间中的守恒律为：

∂U∂t+∂F∂x+∂G∂y=S,\frac{\partial U}{\partial t} + \frac{\partial F}{\partial x} + \frac{\partial G}{\partial y} = S,

经过坐标变换以后，最好仍然保持：

时间积累+计算坐标方向的通量散度=源项\boxed{ \text{时间积累} + \text{计算坐标方向的通量散度} = \text{源项} }

的形式。

通常可以写成类似：

∂(JU)∂τ+∂F^∂ξ+∂G^∂η=JS\boxed{ \frac{\partial(JU)}{\partial \tau} + \frac{\partial\hat F}{\partial\xi} + \frac{\partial\hat G}{\partial\eta} = JS }

其中：

  • JJ 表示几何伸缩；
  • F^,G^\hat F,\hat G 是变换后的通量。

具体表达式会随 Jacobian 定义和是否采用移动网格而有所不同。

这种形式称为强守恒形式。

其意义是：

坐标虽然改变了，但质量、动量和能量仍然应该以真正的通量形式进行离散。

这对有限体积和含激波的可压缩流尤其重要。

静止网格与移动网格

如果网格固定：

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

则计算坐标和物理坐标之间的关系不随时间改变。

但自适应网格中：

x=x(ξ,η,τ),x=x(\xi,\eta,\tau), y=y(ξ,η,τ).y=y(\xi,\eta,\tau).

这里通常令：

τ=t.\tau=t.

需要特别注意：

τ=t\tau=t 并不意味着两个时间偏导算子完全相同。

真正重要的是偏导时保持什么坐标不变。

三种不同的观察方式

固定物理空间

∂∂t∣x,y\left. \frac{\partial}{\partial t} \right|_{x,y}

表示实验室中的固定观察点。

固定计算网格

∂∂τ∣ξ,η\left. \frac{\partial}{\partial\tau} \right|_{\xi,\eta}

表示跟随某个网格节点。

跟随流体质点

DDt\frac{D}{Dt}

表示跟随真实流体粒子。

因此：

物理观察者≠网格观察者≠流体质点观察者\boxed{ \text{物理观察者} \neq \text{网格观察者} \neq \text{流体质点观察者} }

移动网格中的时间变换

若：

x=x(ξ,η,τ),x=x(\xi,\eta,\tau), y=y(ξ,η,τ),y=y(\xi,\eta,\tau),

则固定计算坐标时：

xτ,yτx_\tau,\qquad y_\tau

就是网格在物理空间中的速度。

而控制方程变换时需要：

ξt,ηt,\xi_t,\qquad\eta_t,

它们表示：

固定一个物理点，该点在移动计算坐标系中的坐标变化速度。

两者之间满足：

[ξtηt]=−[xξxηyξyη]−1[xτyτ]\boxed{ \begin{bmatrix} \xi_t\\ \eta_t \end{bmatrix} = - \begin{bmatrix} x_\xi & x_\eta\\ y_\xi & y_\eta \end{bmatrix}^{-1} \begin{bmatrix} x_\tau\\ y_\tau \end{bmatrix} }

因此对于某个流场变量：

ϕ(x,y,t)=Φ(ξ,η,τ),\phi(x,y,t)=\Phi(\xi,\eta,\tau),

有：

∂ϕ∂t∣x,y=∂Φ∂τ∣ξ,η+ξtΦξ+ηtΦη\boxed{ \left. \frac{\partial\phi}{\partial t} \right|_{x,y} = \left. \frac{\partial\Phi}{\partial\tau} \right|_{\xi,\eta} + \xi_t\Phi_\xi + \eta_t\Phi_\eta }

这可以看成“移动坐标系版本的物质导数思想”。

拉伸网格

最简单的非均匀结构网格不是改变网格拓扑，而是改变网格点间距。

其目的通常是：

  • 在壁面附近加密；
  • 在激波附近加密；
  • 在尾迹附近加密；
  • 在远场减少网格点。

基本思想是把有限的计算资源集中到变化最剧烈的区域。

但网格尺寸不能突然变化。

相邻单元尺寸变化过快会：

  • 增大截断误差；
  • 降低迭代稳定性；
  • 使网格质量恶化。

所以拉伸一般要求保持平滑。

贴体网格：椭圆型网格生成

对于任意复杂物体，我们往往只知道边界上的 x,y 而不知道内部网格点应该在哪里。

这正好是一个，已知边界，求内部。

由于椭圆型 PDE 具有全局耦合和平滑性质，可以令坐标函数满足：

∇2ξ=0,∇2η=0.\nabla^2\xi=0, \qquad \nabla^2\eta=0.

或者更加一般的 Poisson 方程。

边界位置由物体几何直接给出，内部坐标则由椭圆方程确定。

为什么选择椭圆型 PDE

椭圆型方程的特点是：整个边界共同影响内部解。

这非常符合网格生成需求。

我们希望：

  • 所有内部网格点同时感受到整个边界形状；
  • 网格线尽可能平滑；
  • 避免局部突然折弯。

因此椭圆网格生成可以理解为：用 PDE 对边界几何进行平滑延拓。

计算形式

实际计算中通常求：

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

原本的 Laplace 方程转换以后可得到类似：

αxξξ−2βxξη+γxηη=0,\alpha x_{\xi\xi} -2\beta x_{\xi\eta} +\gamma x_{\eta\eta}=0, αyξξ−2βyξη+γyηη=0,\alpha y_{\xi\xi} -2\beta y_{\xi\eta} +\gamma y_{\eta\eta}=0,

其中：

α=xη2+yη2,\alpha=x_\eta^2+y_\eta^2, β=xξxη+yξyη,\beta=x_\xi x_\eta+y_\xi y_\eta, γ=xξ2+yξ2.\gamma=x_\xi^2+y_\xi^2.

这里：

β=rξ⋅rη.\beta= \mathbf r_\xi\cdot\mathbf r_\eta.

所以：

β=0\beta=0

意味着两族网格线局部正交。

这些系数不是额外需要求解的未知量，而是由当前网格形状决定的几何量。

自适应网格

固定非均匀网格必须提前猜测哪些区域需要加密。

但流动结构的位置往往会随时间变化。

例如：

  • 激波移动；
  • 涡结构移动；
  • 自由界面移动。

因此可以根据当前流场自动调整网格位置。

典型思想是构造一个 monitor function，例如依赖：

∣∇p∣,∣∇ρ∣,∣ω∣|\nabla p|, \qquad |\nabla\rho|, \qquad |\omega|

等量。

梯度大的地方增加网格密度。

所以自适应网格意味着网格位置也是时间相关未知量。

此时就必须使用的时间度量系数和网格速度（固定观察的计算平面上的一个网格点，在物理平面上运动的速度）。

自适应网格虽然能够节约自由度，但会引入：

  • 网格移动方程；
  • 坐标变换；
  • 网格质量控制；
  • 网格扭曲；
  • 几何守恒问题；
  • 额外通信和数据移动。

分区网格 / 多块结构网格

复杂区域不一定需要用一张网格覆盖。

可以把区域划分成多个：

Ω1,Ω2,…\Omega_1,\Omega_2,\dots

每个区域内部使用自己的结构网格。

这种方法通常称为：

  • zonal grid；
  • multiblock grid。

优点是：

  • 每个局部区域网格容易生成；
  • 可以针对不同区域独立加密；
  • 适合复杂几何；
  • 很适合 domain decomposition。

缺点是：

  • 块间接口必须交换数据；
  • 接口插值可能带来误差；
  • 负载均衡需要额外设计。

从 HPC 角度看，多块网格天然适合一个或若干 block 分配给一个 MPI rank。

因此与并行域分解关系非常紧密。

用于有限体积法的网格

有限差分通常依赖较规则的 stencil，因此传统上更喜欢结构网格。

有限体积法从控制体积分守恒律出发：

ddt∫ΩiU dV+∮∂ΩiF⋅n dS=0.\frac{d}{dt} \int_{\Omega_i}U\,dV + \oint_{\partial\Omega_i} \mathbf F\cdot\mathbf n\,dS =0.

它只需要明确：

  • 每个控制体体积；
  • 每个面的面积和法向；
  • 相邻单元是谁。

因此对网格拓扑的要求更弱。

这使有限体积法非常适合复杂非结构网格

非结构网格

非结构网格可以采用：

二维：

  • 三角形；
  • 四边形；
  • 混合单元。

三维：

  • 四面体；
  • 六面体；
  • 棱柱；
  • 金字塔；
  • 混合网格。

其最大优点是：

可以非常灵活地填充复杂几何。

代价是不能再简单地用：

u[i][j]u[i][j]

描述邻接关系。

通常必须额外保存：

  • cell-to-cell；
  • cell-to-face；
  • face-to-node；
  • boundary type

等拓扑信息。

HPC 视角

结构网格：

(i,j,k)(i,j,k)

访问规则，非常适合：

  • cache；
  • SIMD；
  • stencil；
  • GPU memory coalescing。

非结构网格访问往往需要：

index array\text{index array}

进行间接寻址。

因此：

  • cache locality 更差；
  • SIMD 更困难；
  • GPU memory access 更不规则。

所以几何灵活性和计算规则性之间存在明显 trade-off。

这也是工业 CFD/HPC 中非常重要的问题。

嵌套 / 重叠网格：Chimera Grid

对于多个复杂物体或运动物体，可以分别为每个物体生成一套简单结构网格，然后让不同网格之间相互重叠。

这种方法称为：

overset / Chimera grid\boxed{ \text{overset / Chimera grid} }

基本过程包括：

  1. 每个物体建立自己的局部网格；
  2. 删除落在实体内部的无效单元；
  3. 找到不同网格之间的重叠区域；
  4. 在 donor cell 与 receptor point 之间插值传递流场数据。

优点：

  • 几何处理简单；
  • 运动物体非常方便；
  • 每个局部网格仍可保持高质量结构网格。

缺点：

  • 重叠区域需要插值；
  • 守恒处理比较困难；
  • donor 搜索成本较高；
  • 并行情况下需要大量跨网格通信。

本章方法之间的关系

可以把第 5 章理解成几个逐渐增强的层次：

flowchart TD
    A["规则 Cartesian 网格"] --> B["拉伸结构网格"]
    B --> C["贴体曲线网格"]
    C --> D["自适应移动网格"]
    D --> E["多块 / 分区网格"]

另一条路线则是：

flowchart LR
    A["有限体积"] --> B["非结构网格"]
    B --> C["复杂工业几何"]

以及：

flowchart LR
    A["多个局部结构网格"] --> B["overset / Chimera"]

不同方法主要是在 几何灵活性、网格质量、算法复杂度和计算效率之间做不同取舍。