第4章：欧拉视角（上）

在前两章中，我们从拉格朗日视角探讨了物理模拟，其中计算域随物质运动。本章将转向欧拉视角，在固定的空间网格上求解流体运动方程。我们将深入研究稳定流体算法、压力投影方法、Staggered网格的数学基础，以及高效的线性系统求解器。这些技术构成了现代流体模拟引擎的核心。

4.1 稳定流体与半拉格朗日输送

欧拉方法在固定的空间网格上求解流体方程，避免了拉格朗日方法中的网格变形问题。Stam在1999年提出的"稳定流体"算法通过巧妙结合半拉格朗日输送和隐式粘性求解，实现了无条件稳定的流体模拟，这一突破性工作使得实时流体动画成为可能。

4.1.1 欧拉视角下的Navier-Stokes方程

不可压缩流体的运动由Navier-Stokes方程描述：

$$\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u} + \mathbf{f}$$

$$\nabla \cdot \mathbf{u} = 0$$ 其中：

$\mathbf{u}(\mathbf{x}, t)$：速度场
$p(\mathbf{x}, t)$：压力场
$\rho$：流体密度
$\nu$：运动粘度
$\mathbf{f}$：外力（如重力）

在欧拉框架下，我们在固定点$\mathbf{x}$观察流体性质的变化。物质导数可以展开为： $$\frac{D\phi}{Dt} = \frac{\partial \phi}{\partial t} + (\mathbf{u} \cdot \nabla)\phi$$ 这里$\frac{\partial \phi}{\partial t}$是局部变化率，$(\mathbf{u} \cdot \nabla)\phi$是对流项，描述了物质通过固定点的输送效应。

4.1.2 算子分裂方法

直接求解完整的Navier-Stokes方程非常困难。算子分裂（Operator Splitting）将复杂的PDE分解为多个简单的子问题：

对流步骤： $$\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = 0$$
扩散步骤： $$\frac{\partial \mathbf{u}}{\partial t} = \nu\nabla^2\mathbf{u}$$
外力步骤： $$\frac{\partial \mathbf{u}}{\partial t} = \mathbf{f}$$
投影步骤： $$\frac{\partial \mathbf{u}}{\partial t} = -\frac{1}{\rho}\nabla p, \quad \nabla \cdot \mathbf{u} = 0$$ 每个子步骤的解按顺序组合，得到完整时间步的近似解。这种分裂引入的误差为$O(\Delta t^2)$（Strang分裂）或$O(\Delta t)$（一阶分裂）。

4.1.3 半拉格朗日输送原理

半拉格朗日方法是稳定流体算法的核心创新。考虑对流方程： $$\frac{D\phi}{Dt} = 0$$ 这表示物理量$\phi$沿流体粒子轨迹保持不变。在时刻$t + \Delta t$，位置$\mathbf{x}$处的值等于时刻$t$在上游位置$\mathbf{x}_p$处的值： $$\phi(\mathbf{x}, t + \Delta t) = \phi(\mathbf{x}_p, t)$$ 关键是找到粒子的出发点$\mathbf{x}_p$。一阶近似： $$\mathbf{x}_p = \mathbf{x} - \Delta t \cdot \mathbf{u}(\mathbf{x}, t)$$ 更精确的中点法： $$\mathbf{x}_{mid} = \mathbf{x} - \frac{\Delta t}{2} \cdot \mathbf{u}(\mathbf{x}, t)$$ $$\mathbf{x}_p = \mathbf{x} - \Delta t \cdot \mathbf{u}(\mathbf{x}_{mid}, t)$$ 高阶时间积分方案

Runge-Kutta 4阶（RK4）方法提供更高精度： $$\mathbf{k}_1 = \mathbf{u}(\mathbf{x}, t)$$ $$\mathbf{k}_2 = \mathbf{u}(\mathbf{x} - \frac{\Delta t}{2}\mathbf{k}_1, t + \frac{\Delta t}{2})$$ $$\mathbf{k}_3 = \mathbf{u}(\mathbf{x} - \frac{\Delta t}{2}\mathbf{k}_2, t + \frac{\Delta t}{2})$$ $$\mathbf{k}_4 = \mathbf{u}(\mathbf{x} - \Delta t\mathbf{k}_3, t + \Delta t)$$ $$\mathbf{x}_p = \mathbf{x} - \frac{\Delta t}{6}(\mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3 + \mathbf{k}_4)$$

实践中，2阶中点法通常提供良好的精度-效率平衡。

空间插值方案详解

由于$\mathbf{x}_p$通常不在网格点上，需要插值获得$\phi(\mathbf{x}_p, t)$：

双线性/三线性插值（2D/3D） - 计算简单，$C^0$连续 - 引入一阶数值耗散 - 2D情况：$$\phi(x,y) = (1-s)(1-t)\phi_{00} + s(1-t)\phi_{10} + (1-s)t\phi_{01} + st\phi_{11}$$ 其中$s = (x - x_0)/\Delta x$，$t = (y - y_0)/\Delta y$
Catmull-Rom样条插值 - $C^1$连续，减少数值耗散 - 需要4×4（2D）或4×4×4（3D）网格点 - 可能产生超调（overshoot）
Hermite插值（保持单调性） - 使用WENO限制器防止振荡 - 保持物理量的界限（如密度非负） - 计算成本较高但结果稳定
立方插值配合通量限制器（CIP） - 同时插值函数值和导数 - 紧凑模板，高精度 - 特别适合保持锐利界面

误差分析与改进

半拉格朗日方法的截断误差包含两部分：

时间积分误差：$O(\Delta t^{p+1})$，其中$p$是积分方法的阶数
空间插值误差：$O(\Delta x^{q+1})$，其中$q$是插值的阶数

总误差：$$\epsilon = C_1\Delta t^{p+1} + C_2\frac{|\mathbf{u}|\Delta t}{\Delta x}\Delta x^{q+1}$$ 当$|\mathbf{u}|\Delta t \sim \Delta x$（CFL数约为1）时，空间误差主导。

4.1.4 稳定性分析与CFL条件

传统显式方法受CFL（Courant-Friedrichs-Lewy）条件限制： $$\Delta t \leq \frac{\Delta x}{|\mathbf{u}|_{max}}$$ 这要求时间步长足够小，使得信息在一个时间步内传播不超过一个网格单元。违反CFL条件会导致数值不稳定。

Von Neumann稳定性分析

考虑一维线性对流方程$\frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} = 0$，使用Fourier模式$u = e^{ikx}$分析：

显式迎风格式：放大因子$G = 1 - \sigma(1 - e^{-ik\Delta x})$，其中$\sigma = c\Delta t/\Delta x$ 稳定性要求$|G| \leq 1$，得到$0 \leq \sigma \leq 1$
半拉格朗日格式：放大因子$|G| \leq 1$对所有$\sigma$成立（假设适当的插值）

半拉格朗日方法的革命性优势在于无条件稳定性。即使$\Delta t$很大，方法仍然稳定（尽管精度会降低）。这是因为：

能量耗散：插值引入的数值耗散防止了能量积累
有界性：输送后的值总在输入值的范围内
隐式处理：避免了显式方法的稳定性限制

数值耗散的定量分析

对于不同插值方案，有效耗散系数：

线性插值： $$\nu_{eff} \approx \frac{|\mathbf{u}|\Delta x}{2}\left(1 - \frac{|\mathbf{u}|\Delta t}{\Delta x}\right)\frac{|\mathbf{u}|\Delta t}{\Delta x}$$
立方插值： $$\nu_{eff} \approx \frac{|\mathbf{u}|\Delta x^3}{12\Delta t}\left[\left(\frac{|\mathbf{u}|\Delta t}{\Delta x}\right)^2 - 1\right]^2$$ 当CFL数$\sigma = |\mathbf{u}|\Delta t/\Delta x = 1$时，立方插值的耗散最小。

实际影响与缓解策略

大时间步长($\sigma \gg 1$)会带来：

涡旋衰减：小尺度涡旋能量按$\exp(-2\pi^2\nu_{eff}t/L^2)$衰减
锋面扩散：锐利界面按$\sqrt{\nu_{eff}t}$扩散
相位误差：高频波的传播速度偏差

缓解策略：

自适应时间步长：在需要精度的区域减小$\Delta t$
高阶插值：使用WENO、CIP等方法
反扩散校正：BFECC、MacCormack方法
涡旋限制：检测并增强涡旋结构

高阶输送格式预览

BFECC（Back and Forth Error Compensation and Correction） - 前向-后向输送估计误差 - 二阶精度，减少80%数值耗散 - 计算成本约为标准方法的3倍
MacCormack方法 - 预测-校正思想 - 保持稳定性同时提高精度 - 需要限制器防止新的极值
FLIP（Fluid Implicit Particle）混合 - 结合粒子和网格优势 - 几乎无数值耗散 - 可能引入噪声，需要混合参数调节

4.2 Chorin式压力投影方法

压力投影是不可压缩流体模拟的核心算法。Chorin在1968年提出的投影方法巧妙地将速度场分解为无散度部分和梯度场，通过求解压力泊松方程强制速度场满足不可压缩约束。这一方法的数学优雅性和计算效率使其成为CFD领域的基石。

4.2.1 Helmholtz-Hodge分解

Helmholtz-Hodge定理是投影方法的理论基础：任何光滑向量场$\mathbf{w}$都可以唯一分解为： $$\mathbf{w} = \mathbf{u} + \nabla \phi$$ 其中：

$\mathbf{u}$是无散度场：$\nabla \cdot \mathbf{u} = 0$
$\nabla \phi$是标量势的梯度场

关键性质：

正交性：$\langle \mathbf{u}, \nabla \phi \rangle = 0$（在适当边界条件下）
唯一性：分解在给定边界条件下唯一
投影算子：$\mathbb{P}(\mathbf{w}) = \mathbf{u}$是正交投影

数学上，投影算子可以表示为： $$\mathbb{P} = \mathbb{I} - \nabla(\nabla^2)^{-1}\nabla \cdot$$ 这里$(\nabla^2)^{-1}$表示拉普拉斯算子的逆，需要通过求解泊松方程实现。

4.2.2 压力投影算法推导

从动量方程出发： $$\frac{\partial \mathbf{u}}{\partial t} = -(\mathbf{u} \cdot \nabla)\mathbf{u} - \frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u} + \mathbf{f}$$ 投影方法的核心思想是分两步：

步骤1：计算中间速度场$\mathbf{u}^*$

忽略压力项，求解： $$\frac{\mathbf{u}^* - \mathbf{u}^n}{\Delta t} = -(\mathbf{u}^n \cdot \nabla)\mathbf{u}^n + \nu\nabla^2\mathbf{u}^n + \mathbf{f}$$ 这里$\mathbf{u}^*$通常不满足不可压缩条件：$\nabla \cdot \mathbf{u}^* \neq 0$

步骤2：压力修正

寻找压力$p^{n+1}$使得： $$\frac{\mathbf{u}^{n+1} - \mathbf{u}^*}{\Delta t} = -\frac{1}{\rho}\nabla p^{n+1}$$ 要求$\nabla \cdot \mathbf{u}^{n+1} = 0$，对上式取散度： $$\nabla \cdot \mathbf{u}^{n+1} - \nabla \cdot \mathbf{u}^* = -\frac{\Delta t}{\rho}\nabla^2 p^{n+1}$$ 由于$\nabla \cdot \mathbf{u}^{n+1} = 0$，得到压力泊松方程： $$\nabla^2 p^{n+1} = \frac{\rho}{\Delta t}\nabla \cdot \mathbf{u}^*$$ 求解后，更新速度场： $$\mathbf{u}^{n+1} = \mathbf{u}^* - \frac{\Delta t}{\rho}\nabla p^{n+1}$$

4.2.3 边界条件处理

正确的边界条件对投影方法至关重要：

固体边界（无滑移条件）

速度边界条件：$\mathbf{u} = \mathbf{u}_{wall}$（通常为0）

压力边界条件推导：从动量方程在边界的法向分量： $$\frac{\partial u_n}{\partial t} = -\frac{1}{\rho}\frac{\partial p}{\partial n} + \ldots$$ 由于$u_n$在边界固定，$\frac{\partial u_n}{\partial t} = 0$，得到Neumann条件： $$\frac{\partial p}{\partial n} = \rho \mathbf{n} \cdot \left[(\mathbf{u} \cdot \nabla)\mathbf{u} - \nu\nabla^2\mathbf{u} - \mathbf{f}\right]$$ 实际实现细节：

简化处理：$\frac{\partial p}{\partial n} = 0$（对大多数应用足够精确）

精确处理：需要计算边界上的对流和粘性项：

// 计算边界法向梯度
∂p/∂n = ρ * (u·∇u·n - ν∇²u·n - f·n)

离散化实现：

使用镜像点（ghost cell）方法
一阶精度：$p_{ghost} = p_{interior}$
二阶精度：$p_{ghost} = p_{interior} + \Delta x \cdot \frac{\partial p}{\partial n}$

自由表面边界

压力设为大气压：$p = p_{atm}$（Dirichlet条件）

速度外推：$\frac{\partial \mathbf{u}}{\partial n} = 0$

表面张力处理：考虑表面张力时，压力跳跃： $$p_{liquid} - p_{gas} = \sigma \kappa$$ 其中$\sigma$是表面张力系数，$\kappa$是表面曲率

周期边界 $$\mathbf{u}(x_{min}) = \mathbf{u}(x_{max})$$ $$p(x_{min}) = p(x_{max}) + \Delta p$$ 其中$\Delta p$可能非零（如管道流动）。

实现要点：

修改矩阵结构以包含周期连接
使用循环索引：i_wrap = (i + N) % N
对于压力驱动流，在动量方程中添加$-\nabla p_{drive}$

入口/出口边界

入口（Dirichlet速度）：

指定速度剖面：$\mathbf{u} = \mathbf{u}_{inlet}(y,z,t)$
压力使用Neumann条件：$\frac{\partial p}{\partial n} = 0$

出口（对流出口）：

对流条件：$\frac{\partial \mathbf{u}}{\partial t} + U_{conv}\frac{\partial \mathbf{u}}{\partial n} = 0$
其中$U_{conv}$是对流速度（常取平均流速）
压力：$p = p_{ref}$或$\frac{\partial p}{\partial n} = 0$

浸入边界（Immersed Boundary）

对于复杂几何，使用浸入边界方法：

在笛卡尔网格上定义不规则边界
使用插值或体力方法施加边界条件
适合移动边界问题

数值稳定性考虑：

边界条件的相容性：确保速度和压力边界条件不矛盾例如：不能同时指定速度和压力的Dirichlet条件
角点处理：二维角点、三维棱边需要特殊处理常用平均或优先级方法
时变边界：移动边界需要满足几何守恒律（GCL） $\frac{d}{dt}\int_V dV = \int_S \mathbf{u}_{grid} \cdot \mathbf{n} dS$

4.2.4 压力泊松方程

离散化后的压力泊松方程形成大型稀疏线性系统： $$\mathbf{A}\mathbf{p} = \mathbf{b}$$ 其中：

$\mathbf{A}$是离散拉普拉斯算子
$\mathbf{b} = \frac{\rho}{\Delta t}\nabla \cdot \mathbf{u}^*$

离散拉普拉斯算子详解

2D均匀网格： $$(\nabla^2 p)_{i,j} = \frac{p_{i+1,j} - 2p_{i,j} + p_{i-1,j}}{\Delta x^2} + \frac{p_{i,j+1} - 2p_{i,j} + p_{i,j-1}}{\Delta y^2}$$ 当$\Delta x = \Delta y = h$时，简化为： $$(\nabla^2 p)_{i,j} = \frac{1}{h^2}(p_{i+1,j} + p_{i-1,j} + p_{i,j+1} + p_{i,j-1} - 4p_{i,j})$$
3D情况： $$(\nabla^2 p)_{i,j,k} = \frac{1}{h^2}(p_{i+1,j,k} + p_{i-1,j,k} + p_{i,j+1,k} + p_{i,j-1,k} + p_{i,j,k+1} + p_{i,j,k-1} - 6p_{i,j,k})$$
非均匀网格：使用有限体积方法： $$(\nabla^2 p)_{i,j} = \frac{2}{\Delta x_i(\Delta x_i + \Delta x_{i-1})}p_{i+1,j} + \frac{2}{\Delta x_{i-1}(\Delta x_i + \Delta x_{i-1})}p_{i-1,j} + \ldots$$ 矩阵结构分析

对于$N \times N$网格，系数矩阵$\mathbf{A}$的结构：

[ 4 -1  0 ... -1  0  0 ...]
[-1  4 -1 ...  0 -1  0 ...]
[ 0 -1  4 ...  0  0 -1 ...]
[... ... ... ... ... ... ...]
[-1  0  0 ...  4 -1  0 ...]

特殊性质：

对角占优：$|a_{ii}| \geq \sum_{j \neq i} |a_{ij}|$
对称正定：对于Dirichlet或混合边界条件
条件数：$\kappa \sim O(h^{-2})$，随网格加密恶化
非零元素：2D为5个，3D为7个

兼容性条件与处理

对于纯Neumann边界，泊松方程的解存在当且仅当： $$\int_\Omega \nabla \cdot \mathbf{u}^* dV = \int_{\partial\Omega} \mathbf{u}^* \cdot \mathbf{n} dS = 0$$ 这要求离散速度场满足全局质量守恒。

实际处理方法：

全局质量修正： $$\mathbf{u}^*_{corrected} = \mathbf{u}^* - \frac{1}{|\Omega|}\int_\Omega \nabla \cdot \mathbf{u}^* dV$$ 均匀分布误差，保持局部结构
固定参考点：设$p_{0,0} = 0$，移除常数自由度修改矩阵第一行：$a_{00} = 1, a_{0j} = 0$
投影方法：求解$\mathbf{A}\mathbf{p} = \mathbf{b} - \frac{\mathbf{1}\mathbf{1}^T\mathbf{b}}{N}$ 其中$\mathbf{1}$是全一向量，$N$是未知数个数
自适应兼容性保证：使用保守型离散格式，机器精度内自动满足

效率优化策略

预计算模板：对于固定网格，预计算并存储系数

stencil[0] = 1.0 / (dx * dx)
stencil[1] = 1.0 / (dy * dy)
stencil[2] = -2.0 * (stencil[0] + stencil[1])

边界模板缓存：为不同边界类型预计算修正模板
向量化计算：使用SIMD指令同时处理4/8个网格点

数值精度考虑

舍入误差累积：使用Kahan求和算法计算$\nabla \cdot \mathbf{u}^*$
迭代终止条件：相对残差$|\mathbf{r}|_2 / |\mathbf{b}|_2 < \varepsilon$ 典型值：$\varepsilon = 10^{-6}$到$10^{-8}$
压力单位化：使用$\tilde{p} = p\Delta t/\rho$避免数值问题

4.3 Staggered网格与零空间分析

MAC（Marker-And-Cell）网格是流体模拟的经典数据结构，由Harlow和Welch在1965年提出。通过将速度分量和压力存储在不同位置，Staggered网格自然满足了离散层面的守恒性质，避免了压力场的棋盘模式等数值病态。

4.3.1 MAC网格结构

在MAC网格中，不同物理量存储在网格的不同位置：

2D情况：

压力$p$：存储在网格中心$(i, j)$
$x$方向速度$u$：存储在垂直边中点$(i+\frac{1}{2}, j)$
$y$方向速度$v$：存储在水平边中点$(i, j+\frac{1}{2})$

3D情况：

压力$p$：存储在体元中心$(i, j, k)$
$x$方向速度$u$：存储在$x$垂直面中心$(i+\frac{1}{2}, j, k)$
$y$方向速度$v$：存储在$y$垂直面中心$(i, j+\frac{1}{2}, k)$
$z$方向速度$w$：存储在$z$垂直面中心$(i, j, k+\frac{1}{2})$

这种交错排列的优势：

自然的中心差分：梯度和散度算子无需插值
紧凑的模板：最小的离散化模板
守恒性质：离散算子满足连续算子的关键性质

4.3.2 离散微分算子

散度算子（2D）： $$(\nabla \cdot \mathbf{u})_{i,j} = \frac{u_{i+\frac{1}{2},j} - u_{i-\frac{1}{2},j}}{\Delta x} + \frac{v_{i,j+\frac{1}{2}} - v_{i,j-\frac{1}{2}}}{\Delta y}$$ 注意速度分量正好位于压力单元的边界上，无需插值。

梯度算子（2D）： $$(\nabla p)_x|_{i+\frac{1}{2},j} = \frac{p_{i+1,j} - p_{i,j}}{\Delta x}$$ $$(\nabla p)_y|_{i,j+\frac{1}{2}} = \frac{p_{i,j+1} - p_{i,j}}{\Delta y}$$ 梯度计算位置正好是速度分量的存储位置。

关键性质：离散算子的伴随关系

连续情况下，梯度和负散度互为伴随算子： $$\langle \nabla p, \mathbf{u} \rangle = -\langle p, \nabla \cdot \mathbf{u} \rangle + \text{边界项}$$ MAC网格保持了这一性质（忽略边界）： $$\sum_{i,j} (\nabla p) \cdot \mathbf{u} = -\sum_{i,j} p (\nabla \cdot \mathbf{u})$$ 这保证了能量守恒等重要物理性质。

4.3.3 零空间问题与解决方案

压力零空间

对于纯Neumann边界条件，压力泊松方程： $$\nabla^2 p = f$$ 有一维零空间：常数函数。这是因为$\nabla^2(p + c) = \nabla^2 p$。

解决方案：

固定一点：设$p_{0,0} = 0$
零均值投影：求解后减去平均值
兼容性条件：确保$\int f dV = 0$

速度零空间（棋盘模式）

在同位网格（collocated grid）中，存在棋盘模式：

+1 -1 +1 -1
-1 +1 -1 +1
+1 -1 +1 -1
-1 +1 -1 +1

这种模式下$\nabla \cdot \mathbf{u} = 0$但$\mathbf{u} \neq 0$。

MAC网格通过交错排列自然消除了这种零空间，因为：

散度算子使用相邻速度分量
棋盘模式会产生非零散度

离散Hodge分解的唯一性

在MAC网格上，离散Hodge分解： $$\mathbf{w} = \mathbf{u} + \mathbf{G}p$$ 其中$\mathbf{G}$是离散梯度，$\mathbf{D}$是离散散度，满足：

$\mathbf{D}\mathbf{u} = 0$（无散度）
$\mathbf{D}\mathbf{G} = \mathbf{L}$（离散拉普拉斯）

分解的唯一性（模掉压力常数）由$\text{Ker}(\mathbf{L}) \cap \text{Ker}(\mathbf{G}^T) = \{0\}$保证。

4.3.4 守恒性质分析

MAC网格的一个重要优势是在离散层面保持了关键的守恒定律：

质量守恒

对于不可压缩流体： $$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0$$ 当$\rho$恒定时，简化为$\nabla \cdot \mathbf{u} = 0$。投影步骤精确强制此约束。

动量守恒

考虑无粘性情况，动量方程的积分形式： $$\frac{d}{dt}\int_\Omega \rho \mathbf{u} dV = -\int_{\partial\Omega} p\mathbf{n} dS + \int_\Omega \rho \mathbf{f} dV$$ 离散化保持了这一平衡（数值精度内）。

动能守恒（无粘性）

理想情况下，动能变化率： $$\frac{dE}{dt} = \int_\Omega \mathbf{u} \cdot \mathbf{f} dV$$ 然而，数值方法通常引入耗散。MAC网格配合适当的时间积分可以最小化人工耗散。

涡度守恒

2D无粘性流动中，涡度$\omega = \nabla \times \mathbf{u}$沿流线守恒： $$\frac{D\omega}{Dt} = 0$$ 虽然MAC网格不直接存储涡度，但通过精确的速度插值可以较好保持涡度。

离散分析工具

验证守恒性质的关键是检查离散算子的性质：

反对称性：对流算子应满足$\langle \mathbf{u}, (\mathbf{u} \cdot \nabla)\mathbf{u} \rangle = 0$
正定性：粘性算子应满足$\langle \mathbf{u}, \nabla^2\mathbf{u} \rangle \leq 0$
相容性：离散算子在网格加密时收敛到连续算子

4.4 Krylov子空间求解器

压力泊松方程导出的线性系统通常规模巨大（3D模拟中可达数百万未知数），直接求解器的$O(n^3)$复杂度不可接受。Krylov子空间方法通过在逐步扩大的子空间中寻找近似解，将求解复杂度降至$O(n^{3/2})$或更低，是大规模流体模拟的关键技术。

4.4.1 共轭梯度法原理

共轭梯度法（CG）是求解对称正定线性系统$\mathbf{A}\mathbf{x} = \mathbf{b}$的经典Krylov方法。其核心思想是将求解过程转化为优化问题。

能量泛函最小化

定义能量泛函： $$E(\mathbf{x}) = \frac{1}{2}\mathbf{x}^T\mathbf{A}\mathbf{x} - \mathbf{b}^T\mathbf{x}$$ 其梯度为： $$\nabla E(\mathbf{x}) = \mathbf{A}\mathbf{x} - \mathbf{b} = \mathbf{r}$$ 这里$\mathbf{r}$是残差向量。最小化$E(\mathbf{x})$等价于求解原线性系统。

Krylov子空间

从初始猜测$\mathbf{x}_0$出发，第$k$步的Krylov子空间定义为： $$\mathcal{K}_k(\mathbf{A}, \mathbf{r}_0) = \text{span}\{\mathbf{r}_0, \mathbf{A}\mathbf{r}_0, \mathbf{A}^2\mathbf{r}_0, \ldots, \mathbf{A}^{k-1}\mathbf{r}_0\}$$ 其中$\mathbf{r}_0 = \mathbf{b} - \mathbf{A}\mathbf{x}_0$是初始残差。

共轭方向

CG的关键创新是使用$\mathbf{A}$-共轭（正交）的搜索方向： $$\mathbf{p}_i^T\mathbf{A}\mathbf{p}_j = 0, \quad i \neq j$$ 这保证了：

每个方向只需搜索一次
理论上最多$n$步收敛（$n$为未知数个数）
误差在$\mathbf{A}$-范数下单调递减

算法流程

初始化：
  选择初始猜测 x_0
  r_0 = b - Ax_0
  p_0 = r_0

对于 k = 0, 1, 2, ... 直到收敛：
  α_k = (r_k^T r_k) / (p_k^T A p_k)     # 步长
  x_{k+1} = x_k + α_k p_k               # 更新解
  r_{k+1} = r_k - α_k A p_k             # 更新残差
  β_k = (r_{k+1}^T r_{k+1}) / (r_k^T r_k)  # 共轭系数
  p_{k+1} = r_{k+1} + β_k p_k           # 更新搜索方向

收敛性质

误差的能量范数满足： $$|\mathbf{e}_k|_\mathbf{A} \leq 2\left(\frac{\sqrt{\kappa} - 1}{\sqrt{\kappa} + 1}\right)^k |\mathbf{e}_0|_\mathbf{A}$$ 其中$\kappa = \lambda_{max}/\lambda_{min}$是条件数，$\mathbf{e}_k = \mathbf{x}^* - \mathbf{x}_k$是误差向量。

4.4.2 预条件技术

原始CG的收敛速度强烈依赖于矩阵条件数。预条件通过变换系统降低条件数，显著加速收敛。

预条件的数学原理

给定预条件矩阵$\mathbf{M}$，将原系统变换为： $$\mathbf{M}^{-1}\mathbf{A}\mathbf{x} = \mathbf{M}^{-1}\mathbf{b}$$ 理想情况下，$\mathbf{M} \approx \mathbf{A}$使得$\mathbf{M}^{-1}\mathbf{A} \approx \mathbf{I}$，条件数接近1。

常用预条件器

Jacobi预条件 $$\mathbf{M} = \text{diag}(\mathbf{A})$$

简单高效，易于并行
对角占优系统效果好
存储需求：$O(n)$

不完全Cholesky分解（IC） $$\mathbf{M} = \mathbf{L}\mathbf{L}^T$$ 其中$\mathbf{L}$是稀疏下三角矩阵，保持与$\mathbf{A}$相同的稀疏模式。

效果优于Jacobi
需要符号分解预处理
并行性较差

修正不完全Cholesky（MIC） 通过补偿被丢弃的元素，保持行和： $$\sum_j M_{ij} = \sum_j A_{ij}$$

对泊松方程特别有效
保持了离散拉普拉斯的关键性质

多重网格预条件 使用一个或几个多重网格V循环作为预条件器： $$\mathbf{M}^{-1} \approx \text{MG}(\mathbf{A}, \cdot)$$

最优复杂度$O(n)$
对各种问题鲁棒
实现复杂

预条件CG算法（PCG）

初始化：
  x_0, r_0 = b - Ax_0
  解 Mz_0 = r_0
  p_0 = z_0

迭代：
  α_k = (r_k^T z_k) / (p_k^T A p_k)
  x_{k+1} = x_k + α_k p_k
  r_{k+1} = r_k - α_k A p_k
  解 Mz_{k+1} = r_{k+1}
  β_k = (r_{k+1}^T z_{k+1}) / (r_k^T z_k)
  p_{k+1} = z_{k+1} + β_k p_k

4.4.3 无矩阵方法

在流体模拟中，显式存储系数矩阵$\mathbf{A}$常常不必要且浪费内存。无矩阵方法仅需要计算矩阵-向量乘积$\mathbf{A}\mathbf{v}$。

矩阵-向量乘积的隐式计算

对于泊松方程$\nabla^2 p = f$，矩阵-向量乘积等价于：

对每个网格点 (i,j,k):
  (Ap)[i,j,k] = (6p[i,j,k] - p[i-1,j,k] - p[i+1,j,k]

                          - p[i,j-1,k] - p[i,j+1,k]
                          - p[i,j,k-1] - p[i,j,k+1]) / h²

优势：

内存效率：3D均匀网格仅需存储7个系数（或更少）
缓存友好：避免间接内存访问
易于向量化：规则的计算模式
边界条件灵活：直接在计算中处理

变系数情况

对于变密度流体： $$\nabla \cdot \left(\frac{1}{\rho}\nabla p\right) = \nabla \cdot \mathbf{u}^*$$ 离散化后：

(Ap)[i,j,k] = (p[i+1,j,k] - p[i,j,k]) / (ρ[i+1/2,j,k] h²)

            + (p[i-1,j,k] - p[i,j,k]) / (ρ[i-1/2,j,k] h²)
            + ... (其他方向类似)

性能优化技巧

循环分块：提高缓存局部性
SIMD向量化：使用AVX/SSE指令
GPU实现：高度并行的模板操作
混合精度：单精度计算，双精度累加

4.4.4 收敛性分析

Krylov方法的收敛行为对算法选择和参数调优至关重要。

谱分析

对于对称正定矩阵，CG的收敛速度由特征值分布决定：

聚集的特征值：快速收敛
分散的特征值：缓慢收敛
离群特征值：初期收敛慢，但可通过预条件消除

超线性收敛

CG常表现出超线性收敛：后期收敛速度加快。这是因为：

Ritz值逐渐逼近极端特征值
有效条件数随迭代改善
误差在"困难"的特征方向被消除后，收敛加速

收敛准则

实践中常用的停止准则：

相对残差 $$\frac{|\mathbf{r}_k|_2}{|\mathbf{b}|_2} < \varepsilon$$

直观但可能过于严格
对于病态系统可能永不满足

相对残差下降 $$\frac{|\mathbf{r}_k|_2}{|\mathbf{r}_0|_2} < \varepsilon$$

更实用，但依赖初始猜测

能量范数误差估计 $$\frac{|\mathbf{e}_k|_\mathbf{A}}{|\mathbf{e}_0|_\mathbf{A}} \approx \sqrt{\frac{\mathbf{r}_k^T\mathbf{M}^{-1}\mathbf{r}_k}{\mathbf{r}_0^T\mathbf{M}^{-1}\mathbf{r}_0}}$$

理论上更合理
需要预条件器

其他Krylov方法

BiCGSTAB：适用于非对称系统 - 避免复数运算 - 收敛可能不规则
GMRES：适用于一般非对称系统 - 理论保证收敛 - 内存需求随迭代增长 - 常用重启版本GMRES(m)
MINRES：适用于对称不定系统 - 三项递推关系 - 残差最小化

流体模拟中的选择

不可压缩流：PCG配合MIC或多重网格预条件
可压缩流：GMRES配合ILU预条件
多相流：灵活GMRES(FGMRES)适应变化的预条件器

4.5 多重网格方法

多重网格方法是求解大规模线性系统和偏微分方程的最优算法之一。通过在不同尺度的网格上协同工作，多重网格可以高效消除各种频率的误差分量，实现$O(n)$的最优复杂度。这种方法的数学优雅性和实用性使其成为现代CFD的核心技术。

4.5.1 多重网格基本原理

误差的频谱分析

考虑一维1D泊松方程： $$-\frac{d^2u}{dx^2} = f, \quad x \in [0, 1]$$ 离散化后的误差可以展开为僅量基函数： $$e(x) = \sum_{k=1}^{n} a_k \sin(k\pi x)$$ 不同频率分量的特性：

高频误差（$k \approx n$）：在网格尺度上振荡
低频误差（$k \ll n$）：光滑变化

迭代方法的平滑性质

Jacobi或Gauss-Seidel等经典迭代法对不同频率误差的消除效率：

对于误差分量$e_k = \sin(k\pi x)$，Jacobi迭代的衰减因子： $$\rho_k = 1 - \frac{2\omega}{\lambda_{max}} \cdot 2(1 - \cos(k\pi h))$$ 其中$\omega$是松弛因子，$h = 1/n$是网格宽度。

关键观察：

高频：$|\rho_k| \ll 1$，快速衰减
低频：$|\rho_k| \approx 1$，几乎不衰减

粗网格修正原理

多重网格的核心洞察：细网格上的低频误差在粗网格上变成高频。

考虑波长$\lambda = 2/k$的误差：

细网格（$n$个点）：$\lambda/h = 2n/k$ 个网格点
粗网格（$n/2$个点）：$\lambda/(2h) = n/k$ 个网格点

当$k$较小时，误差在粗网格上相对频率更高，更容易被平滑迭代消除。

两网格校正格式

基本的两网格算法：

预平滑：在细网格上迭代$\nu_1$次 $$\mathbf{x}^h \leftarrow \text{Smooth}^{\nu_1}(\mathbf{A}^h, \mathbf{x}^h, \mathbf{b}^h)$$
残差计算： $$\mathbf{r}^h = \mathbf{b}^h - \mathbf{A}^h\mathbf{x}^h$$
限制：将残差投影到粗网格 $$\mathbf{r}^{2h} = \mathbf{I}_{h}^{2h} \mathbf{r}^h$$
粗网格求解：解误差方程 $$\mathbf{A}^{2h}\mathbf{e}^{2h} = \mathbf{r}^{2h}$$
延拓：将修正插值回细网格 $$\mathbf{x}^h \leftarrow \mathbf{x}^h + \mathbf{I}_{2h}^{h} \mathbf{e}^{2h}$$
后平滑：再次迭代$\nu_2$次 $$\mathbf{x}^h \leftarrow \text{Smooth}^{\nu_2}(\mathbf{A}^h, \mathbf{x}^h, \mathbf{b}^h)$$

4.5.2 限制与延拓算子

限制和延拓算子在不同网格层级之间传递信息，其设计直接影响多重网格的效率和收敛性质。

几何多重网格的标准算子

一维情况

全权重限制（Full weighting）： $$r_{2i}^{2h} = \frac{1}{4}r_{2i-1}^h + \frac{1}{2}r_{2i}^h + \frac{1}{4}r_{2i+1}^h$$ 线性插值延拓： $$e_{2i}^h = e_i^{2h}, \quad e_{2i+1}^h = \frac{1}{2}(e_i^{2h} + e_{i+1}^{2h})$$

二维情况

九点限制（全权重）： $$r_{i,j}^{2h} = \frac{1}{16}\begin{bmatrix} 1 & 2 & 1 \\ 2 & 4 & 2 \\ 1 & 2 & 1 \end{bmatrix} \otimes r^h$$ 双线性插值延拓： $$e_{2i+\alpha,2j+\beta}^h = (1-\alpha)(1-\beta)e_{i,j}^{2h} + \alpha(1-\beta)e_{i+1,j}^{2h} + (1-\alpha)\beta e_{i,j+1}^{2h} + \alpha\beta e_{i+1,j+1}^{2h}$$ 其中$\alpha, \beta \in \{0, 1\}$。

算子的伴随性质

为保持稳定性，限制和延拓算子应满足： $$(\mathbf{I}_{h}^{2h})^T = c \cdot \mathbf{I}_{2h}^{h}$$ 对于上述标准选择，$c = 4$（2D情况）。

Galerkin条件

粗网格算子可通过Galerkin投影自动生成： $$\mathbf{A}^{2h} = \mathbf{I}_{h}^{2h} \mathbf{A}^h \mathbf{I}_{2h}^{h}$$ 这保证了：

粗网格算子继承细网格的对称性
变分性质得到保持
离散化误差最小

选择标准

精度要求：高阶插值减少高频混叠
计算效率：简单算子减少计算量
并行性：局部算子更适合并行
鲁棒性：对各种问题稳定工作

4.5.3 V-循环与W-循环

多重网格的效率很大程度上取决于网格层级的遍历策略。

V-循环

最简单和最常用的策略，形状如字母V：

V-Cycle(A^l, x^l, b^l, l):
  if l == l_max:
    x^l = (A^l)^{-1} b^l  // 直接求解
  else:
    x^l = Smooth^{ν1}(A^l, x^l, b^l)  // 预平滑
    r^l = b^l - A^l x^l
    r^{l+1} = Restrict(r^l)
    e^{l+1} = 0
    V-Cycle(A^{l+1}, e^{l+1}, r^{l+1}, l+1)  // 递归
    x^l = x^l + Prolongate(e^{l+1})
    x^l = Smooth^{ν2}(A^l, x^l, b^l)  // 后平滑

复杂度分析（3D情况）：

第$l$层网格点数：$n/8^l$
总工作量：$O(n)(1 + 1/8 + 1/64 + ...) = O(n)$

W-循环

每个粗网格访问两次，形状如字母W：

W-Cycle(A^l, x^l, b^l, l):
  if l == l_max:
    x^l = (A^l)^{-1} b^l
  else:
    // 与V-循环相同的步骤，但递归调用两次
    W-Cycle(A^{l+1}, e^{l+1}, r^{l+1}, l+1)
    e^{l+1} = e^{l+1} + Δe^{l+1}
    W-Cycle(A^{l+1}, Δe^{l+1}, r^{l+1}, l+1)

F-循环

介于V和W之间的折中方案：

第一次下降使用V-循环
第二次下降使用W-循环

收敛性能比较

循环类型	每次迭代工作量	收敛因子	适用场景
V-循环	$O(n)$	0.1-0.2	标准选择，效率高
W-循环	$O(n\log n)$	0.05-0.1	困难问题，更鲁棒
F-循环	$O(n)$	0.08-0.15	平衡效率和鲁棒性

Full Multigrid (FMG)

使用粗网格提供良好初始猜测：

在最粗网格求解
插值到细一级网格
执行V-循环改进
重复直到最细网格

FMG可以在$O(n)$工作量内达到离散化误差级别的精度。

4.5.4 代数多重网格（AMG）

几何多重网格需要规则网格结构，而代数多重网格直接基于矩阵的代数性质构造多层级结构，适用于非结构网格和复杂问题。

AMG的核心思想

“光滑误差”的代数特征： $$\mathbf{A}\mathbf{e} \approx 0$$ 这意味着误差在代数意义上“光滑”。AMG通过分析矩阵元素间的强耦合关系构建粗网格。

强耦合的定义

变量$i$和$j$之间的强耦合条件： $$|a_{ij}| \geq \theta \max_{k \neq i} |a_{ik}|$$ 典型阈值$\theta = 0.25$。强耦合意味着两个变量紧密相关。

粗网格选择算法

第一阶段：分类变量 - C-点（粗网格点）：保留到粗网格 - F-点（细网格点）：仅在细网格存在
选择策略： - 每个F-点至少强耦合于一个C-点 - C-点之间避免强耦合 - 最大化C-点数量（保持粗网格有效性）

插值算子构造

代数插值基于最小化能量的原则：

对于F-点$i$，其值由强耦合的C-点插值： $$e_i = \sum_{j \in C_i} w_{ij} e_j$$ 权重计算（经典AMG）： $$w_{ij} = -\frac{a_{ij} + \sum_{k \in F_i^s} \frac{a_{ik}a_{kj}}{\sum_{m \in C_i} a_{km}}}{a_{ii} + \sum_{k \in F_i^w} a_{ik}}$$ 其中：

$C_i$：与$i$强耦合的C-点
$F_i^s$：与$i$强耦合的F-点
$F_i^w$：与$i$弱耦合的F-点

AMG的优势与挑战

优势：

通用性：不依赖网格结构
自适应：自动适应问题特性
鲁棒性：对各类问题有效
可扩展性：大规模并行性能好

挑战：

设置成本：构建多层级结构耗时
内存需求：存储插值算子
参数敏感：需要调优阈值等参数
理论保证：收敛性分析复杂

平滑者的选择

AMG中常用的平滑者：

Jacobi：并行性好，但效果一般
Gauss-Seidel：效果好，但串行
多色GS：平衡并行性和效果
ILU平滑：强大但昂贵

流体模拟中的应用

压力泊松方程的AMG求解：

对非结构网格特别有效
处理变系数和各向异性
适应复杂边界条件
与预条件Krylov方法结合使用

本章小结

本章深入探讨了欧拉视角下的流体模拟核心技术。我们从稳定流体算法开始，逐步深入到数值方法的数学基础和工程实现。以下是本章的关键要点：

核心概念

半拉格朗日输送 - 无条件稳定性是其最大优势 - 通过回溯粒子轨迹实现对流项的隐式处理 - 数值耗散是代价，需要高阶方法改进
压力投影方法 - Helmholtz-Hodge分解提供理论基础 - 两步法：先计算中间速度，再通过压力修正 - 泊松方程是计算瓶颈
Staggered网格 - 消除棋盘模式等数值病态 - 保持离散层面的守恒性质 - 中心差分无需插值
高效线性求解器 - Krylov子空间方法：$O(n^{3/2})$复杂度 - 预条件技术显著加速收敛 - 多重网格方法：最优$O(n)$复杂度

关键公式总结

Navier-Stokes方程 $$\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u} + \mathbf{f}$$ $$\nabla \cdot \mathbf{u} = 0$$
半拉格朗日输送 $$\phi(\mathbf{x}, t + \Delta t) = \phi(\mathbf{x}_p, t), \quad \mathbf{x}_p = \mathbf{x} - \Delta t \cdot \mathbf{u}(\mathbf{x}, t)$$
压力泊松方程 $$\nabla^2 p^{n+1} = \frac{\rho}{\Delta t}\nabla \cdot \mathbf{u}^*$$
CG收敛率 $$|\mathbf{e}_k|_\mathbf{A} \leq 2\left(\frac{\sqrt{\kappa} - 1}{\sqrt{\kappa} + 1}\right)^k |\mathbf{e}_0|_\mathbf{A}$$
多重网格限制/延拓 $$\mathbf{r}^{2h} = \mathbf{I}_{h}^{2h} \mathbf{r}^h, \quad \mathbf{x}^h \leftarrow \mathbf{x}^h + \mathbf{I}_{2h}^{h} \mathbf{e}^{2h}$$

实践要点

算法选择 - 实时应用：稳定流体 + Jacobi预条件CG - 高精度模拟：BFECC/MacCormack + 多重网格 - 复杂几何：非结构网格 + AMG
性能优化 - 无矩阵方法减少内存带宽 - SIMD向量化加速模板计算 - 领域分解实现大规模并行
数值稳定性 - CFL条件虽不是必须，但影响精度 - 预条件器的选择对收敛速度至关重要 - 边界条件处理需要特别注意兼容性

与拉格朗日方法的对比

| 特性 | 欧拉方法 | 拉格朗日方法 |

特性	欧拉方法	拉格朗日方法
网格结构	固定，规则	随物质运动
边界处理	简单	复杂
大变形	无问题	可能失效
数值耗散	较大	较小
计算效率	高	中等
适用范围	流体、烟雾	固体、布料

欧拉方法特别适合大规模流体模拟，而下一章我们将继续探讨高级输送格式和自由表面追踪等进阶技术。