第六章：高级输送格式与等势面方法

在前面的章节中，我们学习了欧拉视角下的基本输送格式和简单的自由表面处理方法。然而，基础的半拉格朗日方法存在较大的数值耗散，简单的标记方法难以精确描述复杂的界面演化。本章将深入探讨高阶输送格式，介绍基于有符号距离场的等势面方法，以及相关的渲染技术。这些高级技术是构建高质量流体仿真系统的关键组件。

本章的学习目标包括：

• 掌握WENO等高阶输送格式的原理和实现
• 理解有符号距离场(SDF)的构建和维护方法
• 熟练运用Level Set方法进行界面追踪
• 了解各种自由表面追踪技术的优缺点
• 掌握基于物理仿真数据的渲染技术

6.1 高阶输送格式深入

在第四章中，我们介绍了基本的半拉格朗日输送方法。虽然该方法无条件稳定，但存在严重的数值耗散问题。为了在保持稳定性的同时减少数值耗散，研究者们开发了一系列高阶输送格式。

6.1.1 WENO格式

WENO (Weighted Essentially Non-Oscillatory) 格式是一类高精度的数值格式，特别适合处理含有激波或不连续的问题。其核心思想是通过多个子模板的加权组合来重构数值通量，权重根据局部光滑度动态调整。

对于标量输送方程 $\frac{\partial u}{\partial t} + \frac{\partial f(u)}{\partial x} = 0$，五阶WENO格式使用三个三点子模板：

$$S_0 = \{x_{i-2}, x_{i-1}, x_i\}, \quad S_1 = \{x_{i-1}, x_i, x_{i+1}\}, \quad S_2 = \{x_i, x_{i+1}, x_{i+2}\}$$ 每个子模板上的三阶多项式重构为： $$p_k(x) = u_i + \frac{u'_i}{1!}(x-x_i) + \frac{u''_i}{2!}(x-x_i)^2$$ 光滑度指标衡量每个子模板上解的光滑程度： $$\beta_0 = \frac{13}{12}(u_{i-2} - 2u_{i-1} + u_i)^2 + \frac{1}{4}(u_{i-2} - 4u_{i-1} + 3u_i)^2$$

$$\beta_1 = \frac{13}{12}(u_{i-1} - 2u_i + u_{i+1})^2 + \frac{1}{4}(u_{i-1} - u_{i+1})^2$$

$$\beta_2 = \frac{13}{12}(u_i - 2u_{i+1} + u_{i+2})^2 + \frac{1}{4}(3u_i - 4u_{i+1} + u_{i+2})^2$$ 非线性权重通过以下公式计算： $$\alpha_k = \frac{d_k}{(\epsilon + \beta_k)^2}, \quad \omega_k = \frac{\alpha_k}{\sum_{j=0}^2 \alpha_j}$$ 其中 $d_0 = 1/10$, $d_1 = 6/10$, $d_2 = 3/10$ 是理想权重，$\epsilon = 10^{-6}$ 防止除零。

最终的数值通量为： $$\hat{f}_{i+1/2} = \sum_{k=0}^2 \omega_k f^{(k)}_{i+1/2}$$

6.1.2 TVD限制器

总变差递减 (Total Variation Diminishing, TVD) 条件是保证数值格式稳定性的重要准则。对于标量守恒律，TVD条件要求： $$TV(u^{n+1}) \leq TV(u^n)$$ 其中总变差定义为： $$TV(u) = \sum_i |u_{i+1} - u_i|$$ TVD条件保证了数值解不会产生非物理的振荡。为了构造TVD格式，我们使用通量限制器。考虑通量形式： $$f_{i+1/2} = f^L_{i+1/2} + \phi(r_i)(f^H_{i+1/2} - f^L_{i+1/2})$$ 其中 $f^L$ 是低阶通量（如迎风格式），$f^H$ 是高阶通量（如中心差分），$\phi(r)$ 是限制器函数，$r$ 是连续性指标： $$r_i = \frac{u_i - u_{i-1}}{u_{i+1} - u_i}$$

6.1.3 通量限制器

常用的通量限制器包括：

minmod限制器： $$\phi_{minmod}(r) = \max(0, \min(1, r))$$ 这是最保守的限制器，提供强单调性保证但可能过度耗散。

superbee限制器： $$\phi_{superbee}(r) = \max(0, \min(1, 2r), \min(2, r))$$ superbee限制器较为激进，能更好地保持间断但可能在光滑区域产生振荡。

van Leer限制器： $$\phi_{vanLeer}(r) = \frac{r + |r|}{1 + |r|}$$ van Leer限制器提供了良好的平衡，在光滑区域达到二阶精度。

MC (Monotonized Central)限制器： $$\phi_{MC}(r) = \max(0, \min((1+r)/2, 2, 2r))$$ 所有限制器都必须满足以下条件以保证TVD性质：

• $\phi(r) = 0$ 当 $r < 0$（局部极值处）
• $0 \leq \phi(r) \leq 2$
• $0 \leq \phi(r)/r \leq 2$（当 $r > 0$）

6.1.4 特征线方法

特征线方法基于这样的观察：沿着特征线，偏微分方程可以化为常微分方程。对于输送方程： $$\frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} = 0$$ 特征线定义为： $$\frac{dx}{dt} = c$$ 沿特征线，解保持常数： $$\frac{du}{dt} = \frac{\partial u}{\partial t} + \frac{dx}{dt}\frac{\partial u}{\partial x} = \frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} = 0$$ 对于非线性问题 $\frac{\partial u}{\partial t} + \frac{\partial f(u)}{\partial x} = 0$，特征速度为： $$\lambda = \frac{df}{du}$$ 在数值实现中，我们可以：

1. 追踪特征线：从网格点 $(x_i, t^{n+1})$ 向后追踪到时间 $t^n$
2. 插值求值：在追踪点处插值得到解的值
3. 更新解：将插值得到的值赋给新时刻的解

这种方法的优势在于能够精确处理线性输送，并且可以使用高阶插值来提高精度。对于非线性问题，可以使用局部线性化或Riemann求解器。

6.2 有符号距离场(SDF)

有符号距离场是描述隐式曲面的强大工具，在计算机图形学和计算流体力学中有广泛应用。SDF不仅能精确表示复杂的几何形状，还提供了丰富的几何信息，如法向量和曲率。

6.2.1 距离场的定义与性质

有符号距离函数 $\phi(\mathbf{x})$ 定义为点 $\mathbf{x}$ 到最近界面的有符号距离： $$\phi(\mathbf{x}) = \begin{cases} +d(\mathbf{x}, \partial\Omega) & \text{如果 } \mathbf{x} \in \Omega^+ \text{（外部）} \\ 0 & \text{如果 } \mathbf{x} \in \partial\Omega \text{（界面上）} \\ -d(\mathbf{x}, \partial\Omega) & \text{如果 } \mathbf{x} \in \Omega^- \text{（内部）} \end{cases}$$ 其中 $d(\mathbf{x}, \partial\Omega) = \min_{\mathbf{y} \in \partial\Omega} ||\mathbf{x} - \mathbf{y}||$ 是到界面的欧几里得距离。

SDF具有以下重要性质：

1. 梯度性质： $$||\nabla\phi|| = 1 \quad \text{（几乎处处成立）}$$ 这个性质称为Eikonal方程，表明SDF的梯度模长为1。证明：考虑两个相邻点 $\mathbf{x}_1$ 和 $\mathbf{x}_2$，设它们在界面上的最近点分别为 $\mathbf{p}_1$ 和 $\mathbf{p}_2$。根据三角不等式： $$|\phi(\mathbf{x}_1) - \phi(\mathbf{x}_2)| \leq ||\mathbf{x}_1 - \mathbf{x}_2||$$ 因此 $||\nabla\phi|| \leq 1$。在非奇异点处，等号成立。

2. 法向量计算： $$\mathbf{n} = \frac{\nabla\phi}{||\nabla\phi||} = \nabla\phi$$ 由于 $||\nabla\phi|| = 1$，梯度直接给出了单位外法向量。

3. 曲率计算： 平均曲率可以通过以下公式计算： $$\kappa = \nabla \cdot \mathbf{n} = \nabla \cdot \nabla\phi = \Delta\phi$$ 对于二维情况： $$\kappa = \frac{\phi_{xx}\phi_y^2 - 2\phi_x\phi_y\phi_{xy} + \phi_{yy}\phi_x^2}{(\phi_x^2 + \phi_y^2)^{3/2}}$$ 当 $||\nabla\phi|| = 1$ 时，简化为： $$\kappa = \phi_{xx} + \phi_{yy}$$

6.2.2 快速行进法(FMM)

快速行进法是求解Eikonal方程的高效算法，时间复杂度为 $O(N \log N)$。考虑更一般的Eikonal方程： $$||\nabla T(\mathbf{x})|| = \frac{1}{F(\mathbf{x})}, \quad T|_{\Gamma} = 0$$ 其中 $F(\mathbf{x}) > 0$ 是速度函数，$\Gamma$ 是初始界面。当 $F = 1$ 时，$T$ 就是距离函数。

FMM算法基于Dijkstra最短路径算法的思想，维护三类网格点：

• Known：已确定最终值的点
• Trial：邻接Known点，值可能更新的点
• Far：尚未处理的点

算法步骤：

1. 初始化： - 将界面点标记为Known，设置 $T = 0$ - 将界面的邻居点标记为Trial，计算初始值 - 其余点标记为Far，设置 $T = \infty$

2. 主循环：

while Trial集合非空:
    选择Trial中T值最小的点，移入Known
    更新该点的所有Far邻居，移入Trial
    更新该点的所有Trial邻居的T值

3. 局部更新公式（二维情况）： 对于点 $(i,j)$，使用迎风差分离散Eikonal方程： $$\max(D^{-x}_{ij}T, -D^{+x}_{ij}T, 0)^2 + \max(D^{-y}_{ij}T, -D^{+y}_{ij}T, 0)^2 = \frac{1}{F_{ij}^2}$$ 这导出二次方程： $$aT_{ij}^2 - 2bT_{ij} + c = 0$$ 其中：

• $a = n$ （使用的邻居数）
• $b = \sum T_{neighbor}$
• $c = \sum T_{neighbor}^2 - h^2/F_{ij}^2$

6.2.3 快速扫描法(FSM)

快速扫描法是另一种求解Eikonal方程的方法，使用Gauss-Seidel迭代和特定的扫描顺序。虽然理论复杂度为 $O(N)$，但常数因子较大。

FSM的核心思想是利用特征线的传播方向。在2D情况下，使用4个扫描方向：

1. 从左下到右上：$i = 1:I, j = 1:J$
2. 从右下到左上：$i = I:1, j = 1:J$
3. 从右上到左下：$i = I:1, j = J:1$
4. 从左上到右下：$i = 1:I, j = J:1$

在3D情况下需要8个扫描方向。每次扫描使用相同的更新公式：

for 每个扫描方向:
    for 按该方向遍历所有点(i,j):
        使用邻居值更新T[i,j]

FSM的优势：

• 实现简单，无需复杂的数据结构
• 内存访问模式规则，缓存友好
• 易于并行化

FSM的劣势：

• 需要多次扫描才能收敛
• 对于复杂几何可能需要更多迭代

6.2.4 重新初始化

在Level Set方法的演化过程中，距离函数性质 $||\nabla\phi|| = 1$ 会逐渐丧失。为了维持数值稳定性和精度，需要定期重新初始化。

重新初始化方程： $$\frac{\partial\phi}{\partial\tau} + S(\phi_0)(||\nabla\phi|| - 1) = 0$$ 其中 $\tau$ 是虚拟时间，$S(\phi_0)$ 是光滑的符号函数： $$S(\phi_0) = \frac{\phi_0}{\sqrt{\phi_0^2 + (\Delta x)^2}}$$ 或者使用更精确的形式： $$S(\phi_0) = \begin{cases} -1 & \text{if } \phi_0 < -\epsilon \\ \phi_0/\epsilon + \phi_0^3/\epsilon^3 & \text{if } |\phi_0| \leq \epsilon \\ 1 & \text{if } \phi_0 > \epsilon \end{cases}$$ 数值实现时，使用迎风格式： $$\frac{\partial\phi}{\partial\tau} = -S(\phi_0)G$$ 其中： $$G = \begin{cases} \sqrt{\max(a^2, b^2) + \max(c^2, d^2)} - 1 & \text{if } S(\phi_0) > 0 \\ \sqrt{\max(A^2, B^2) + \max(C^2, D^2)} - 1 & \text{if } S(\phi_0) < 0 \end{cases}$$ 这里：

• $a = \max(D^{-x}\phi, 0)$, $b = \min(D^{+x}\phi, 0)$
• $c = \max(D^{-y}\phi, 0)$, $d = \min(D^{+y}\phi, 0)$
• $A = \min(D^{-x}\phi, 0)$, $B = \max(D^{+x}\phi, 0)$
• $C = \min(D^{-y}\phi, 0)$, $D = \max(D^{+y}\phi, 0)$

重新初始化的停止准则：

• 固定迭代次数（如5-10次）
• 残差小于阈值：$\max ||\nabla\phi|| - 1| < \epsilon$
• 界面移动量小于阈值：$\max |\phi^{new} - \phi^{old}|_{at\ \phi=0} < \epsilon$

6.3 等势面(Level Set)方法

Level Set方法是追踪移动界面的强大工具，由Osher和Sethian在1988年提出。该方法将界面隐式地表示为高一维函数的零等势面，自然处理拓扑变化，避免了显式界面追踪的复杂性。

6.3.1 界面演化方程

考虑随时间演化的界面 $\Gamma(t)$，我们将其表示为函数 $\phi(\mathbf{x}, t)$ 的零等势面： $$\Gamma(t) = \{\mathbf{x} : \phi(\mathbf{x}, t) = 0\}$$ 对于界面上的任意点 $\mathbf{x}(t)$，有： $$\phi(\mathbf{x}(t), t) = 0$$ 对时间求导得到： $$\frac{\partial\phi}{\partial t} + \nabla\phi \cdot \frac{d\mathbf{x}}{dt} = 0$$ 设界面的法向速度为 $V_n$（向外为正），则： $$\frac{d\mathbf{x}}{dt} \cdot \mathbf{n} = V_n$$ 其中 $\mathbf{n} = \nabla\phi / ||\nabla\phi||$ 是单位外法向量。因此： $$\frac{\partial\phi}{\partial t} + V_n ||\nabla\phi|| = 0$$ 这就是Level Set演化的基本方程。对于更一般的速度场 $\mathbf{V}$： $$\frac{\partial\phi}{\partial t} + \mathbf{V} \cdot \nabla\phi = 0$$ 这是一个Hamilton-Jacobi型方程，需要使用特殊的数值格式求解。

数值格式：

对于输送项 $\mathbf{V} \cdot \nabla\phi$，使用迎风格式： $$(\mathbf{V} \cdot \nabla\phi)_{ij} = \max(u_{ij}, 0)D^{-x}\phi + \min(u_{ij}, 0)D^{+x}\phi + \max(v_{ij}, 0)D^{-y}\phi + \min(v_{ij}, 0)D^{+y}\phi$$ 对于Hamilton-Jacobi项 $V_n||\nabla\phi||$，使用Godunov格式：

当 $V_n > 0$（界面向外扩张）： $$||\nabla\phi|| \approx \sqrt{\max(D^{-x}\phi, 0)^2 + \min(D^{+x}\phi, 0)^2 + \max(D^{-y}\phi, 0)^2 + \min(D^{+y}\phi, 0)^2}$$ 当 $V_n < 0$（界面向内收缩）： $$||\nabla\phi|| \approx \sqrt{\min(D^{-x}\phi, 0)^2 + \max(D^{+x}\phi, 0)^2 + \min(D^{-y}\phi, 0)^2 + \max(D^{+y}\phi, 0)^2}$$

6.3.2 速度延拓

在许多应用中，速度只在界面附近有定义。为了稳定地演化Level Set函数，需要将速度延拓到整个计算域。速度延拓应满足：

1. 在界面上保持原速度值
2. 沿法向保持常数：$\nabla\phi \cdot \nabla V_{ext} = 0$

这确保了延拓速度不会在法向产生额外的变化。延拓方程为： $$\frac{\partial V_{ext}}{\partial\tau} + S(\phi)\mathbf{n} \cdot \nabla V_{ext} = 0$$ 其中 $\tau$ 是虚拟时间，$S(\phi)$ 是符号函数。数值实现时： $$\frac{\partial V_{ext}}{\partial\tau} = -S(\phi)[\max(\mathbf{n}_x, 0)D^{-x}V + \min(\mathbf{n}_x, 0)D^{+x}V + \max(\mathbf{n}_y, 0)D^{-y}V + \min(\mathbf{n}_y, 0)D^{+y}V]$$

6.3.3 质量守恒问题

Level Set方法的主要缺点是不能精确保持质量守恒。这是因为：

1. 数值耗散：数值格式引入的人工粘性导致界面变得模糊
2. 重新初始化误差：重新初始化过程可能移动界面位置
3. 欠采样：网格分辨率不足时，小特征可能消失

质量损失的定量分析： 设初始体积为 $V_0 = \int_{\phi<0} d\mathbf{x}$，经过时间 $t$ 后的体积为 $V(t)$。质量损失率为： $$\frac{dV}{dt} = -\int_{\phi=0} V_n dS$$ 数值误差导致的额外质量损失约为 $O(\Delta x)$。

改进方法：

1. 高阶格式：使用WENO、ENO等高阶格式减少数值耗散
2. 自适应网格：在界面附近加密网格
3. 混合方法：结合粒子方法（如PLS）或VOF方法

6.3.4 粒子等势面(PLS)

粒子等势面方法通过在界面附近撒布拉格朗日粒子来改善质量守恒。基本思想是：

1. 粒子初始化：在界面两侧的窄带内随机撒布粒子 - 正粒子（$\phi > 0$ 区域） - 负粒子（$\phi < 0$ 区域）

2. 粒子输送：使用速度场输送粒子 $$\frac{d\mathbf{x}_p}{dt} = \mathbf{V}(\mathbf{x}_p, t)$$

3. 误差检测：检查粒子是否逃逸到错误的一侧 - 正粒子进入 $\phi < 0$ 区域 → 界面应该向外扩张 - 负粒子进入 $\phi > 0$ 区域 → 界面应该向内收缩

4. Level Set修正：根据逃逸粒子修正 $\phi$ 值

对于每个逃逸粒子 $p$，定义球形修正函数： $$\phi_p(\mathbf{x}) = s_p(r_p - ||\mathbf{x} - \mathbf{x}_p||)$$ 其中 $s_p = \pm 1$ 是粒子符号，$r_p$ 是粒子半径（通常为 $0.5\Delta x$）。

最终的修正值： $$\phi^{corrected} = \begin{cases} \max(\phi, \max_p \phi_p) & \text{正粒子修正} \\ \min(\phi, \min_p \phi_p) & \text{负粒子修正} \end{cases}$$

5. 粒子重采样：定期删除远离界面的粒子，在界面附近添加新粒子

PLS方法的优势：

• 显著改善质量守恒
• 保持小尺度特征
• 实现相对简单

PLS方法的劣势：

• 增加计算和存储开销
• 粒子分布可能不均匀
• 需要处理粒子的添加和删除

6.4 自由表面追踪

自由表面是流体与空气（或真空）的界面，其追踪是流体仿真中的核心问题。除了Level Set方法，还有其他几种重要的界面追踪技术，每种方法都有其独特的优势和局限性。

6.4.1 VOF方法对比

Volume of Fluid (VOF) 方法使用体积分数 $f$ 来表示每个网格单元中流体的占比： $$f_{ij} = \frac{V_{fluid}}{V_{cell}}$$ 其中 $f = 1$ 表示完全充满流体，$f = 0$ 表示完全是空气，$0 < f < 1$ 表示包含界面。

VOF vs Level Set对比：

| 特性 | VOF方法 | Level Set方法 |

VOF的演化方程： $$\frac{\partial f}{\partial t} + \nabla \cdot (f\mathbf{V}) = 0$$ 数值求解的关键挑战是保持 $f$ 的锐利性（避免数值扩散）和有界性（$0 \leq f \leq 1$）。

6.4.2 界面重构技术

VOF方法需要从体积分数重构出界面的几何形状。常用的重构方法包括：

SLIC (Simple Line Interface Calculation)： 界面用与坐标轴平行的直线表示。对于2D情况，界面方向基于邻近单元的 $f$ 值梯度： $$\text{如果} |f_{i+1,j} - f_{i-1,j}| > |f_{i,j+1} - f_{i,j-1}|, \text{则界面垂直}$$ 界面位置通过体积守恒确定： $$x_{interface} = x_i - \frac{\Delta x}{2} + f_{ij} \Delta x$$ PLIC (Piecewise Linear Interface Calculation)： 界面用任意方向的直线段表示。界面法向量通过梯度估计： $$\mathbf{n} = -\frac{\nabla f}{||\nabla f||}$$ 界面方程为 $\mathbf{n} \cdot (\mathbf{x} - \mathbf{x}_0) = 0$，其中 $\mathbf{x}_0$ 通过体积约束确定： $$\int_{cell \cap \{\mathbf{n} \cdot (\mathbf{x} - \mathbf{x}_0) < 0\}} d\mathbf{x} = f_{ij} \cdot V_{cell}$$ 高阶重构： 使用二次或三次多项式重构界面，提高精度但增加复杂度。例如，抛物线重构： $$y = ax^2 + bx + c$$ 系数通过最小二乘拟合邻近单元的界面位置确定。

6.4.3 表面张力计算

表面张力是维持液滴、气泡等形状的重要力。在连续介质框架下，表面张力表现为压力跳变： $$[p] = \sigma \kappa$$ 其中 $\sigma$ 是表面张力系数，$\kappa$ 是界面曲率，$[p]$ 表示跨界面的压力跳变。

CSF (Continuum Surface Force) 模型： 将表面张力转换为体积力： $$\mathbf{F}_{st} = \sigma \kappa \delta(\phi) \nabla \phi$$ 其中 $\delta(\phi)$ 是Dirac delta函数的光滑近似： $$\delta_\epsilon(\phi) = \begin{cases} \frac{1}{2\epsilon}(1 + \cos(\pi\phi/\epsilon)) & |\phi| < \epsilon \\ 0 & |\phi| \geq \epsilon \end{cases}$$ 曲率计算（二维）： $$\kappa = \nabla \cdot \mathbf{n} = \nabla \cdot \left(\frac{\nabla\phi}{||\nabla\phi||}\right)$$ 数值实现： $$\kappa_{ij} = \frac{1}{||\nabla\phi||_{ij}} \left[\frac{\phi_{i+1,j} - \phi_{i-1,j}}{2\Delta x} \cdot \frac{n^x_{i+1,j} - n^x_{i-1,j}}{2\Delta x} + \frac{\phi_{i,j+1} - \phi_{i,j-1}}{2\Delta y} \cdot \frac{n^y_{i,j+1} - n^y_{i,j-1}}{2\Delta y}\right]$$ Ghost Fluid方法中的表面张力： 在压力投影步骤中直接施加压力跳变条件： $$p^+ - p^- = \sigma \kappa$$ 这避免了将表面力转换为体积力带来的数值抹平。

6.4.4 拓扑变化处理

拓扑变化（如液滴合并、断裂）是自由表面流动的重要特征。

Level Set的自然处理： Level Set方法自动处理拓扑变化，无需特殊算法。当两个界面靠近时： $$\phi_{merged} = \min(\phi_1, \phi_2)$$ 对于分离，当界面变薄时自然断开。

VOF的拓扑处理： VOF方法需要显式检测和处理拓扑变化：

1. 合并检测：当两个充满流体的单元相邻时
2. 分离检测：当连通域分裂时（需要连通性分析）
3. 薄膜破裂：当界面厚度小于阈值时强制断开

混合方法CLSVOF： 结合Level Set和VOF的优势：

• 使用Level Set计算几何信息
• 使用VOF保证质量守恒
• 通过VOF约束修正Level Set： $$\int_{cell} H(-\phi) d\mathbf{x} = f_{ij} \cdot V_{cell}$$ 其中 $H$ 是Heaviside函数。

6.5 渲染技术

物理仿真的最终目的往往是生成视觉效果。本节介绍几种重要的渲染技术，特别是如何高效地渲染流体表面和体积数据。

6.5.1 路径追踪(Path Tracing)

路径追踪是基于物理的渲染方法，通过模拟光线在场景中的传播来生成逼真的图像。渲染方程（Kajiya 1986）描述了光的传输： $$L_o(\mathbf{x}, \omega_o) = L_e(\mathbf{x}, \omega_o) + \int_{\Omega} f_r(\mathbf{x}, \omega_i, \omega_o) L_i(\mathbf{x}, \omega_i) (\omega_i \cdot \mathbf{n}) d\omega_i$$ 其中：

• $L_o$：出射辐射度
• $L_e$：自发光
• $f_r$：双向反射分布函数(BRDF)
• $L_i$：入射辐射度
• $\Omega$：半球立体角

Monte Carlo积分： 由于渲染方程的积分难以解析求解，使用Monte Carlo方法： $$L_o \approx L_e + \frac{1}{N} \sum_{i=1}^N \frac{f_r(\omega_i) L_i(\omega_i) (\omega_i \cdot \mathbf{n})}{p(\omega_i)}$$ 其中 $p(\omega_i)$ 是采样概率密度函数。

俄罗斯轮盘赌： 为了保证算法收敛，使用俄罗斯轮盘赌终止光线：

if (random() < p_rr):
    继续追踪，权重乘以 1/p_rr
else:
    终止光线

重要性采样： 根据BRDF的形状采样方向，提高收敛速度。对于漫反射表面，使用余弦加权采样： $$p(\theta, \phi) = \frac{\cos\theta}{\pi}$$ 生成采样方向： $$\theta = \arccos(\sqrt{1 - \xi_1}), \quad \phi = 2\pi\xi_2$$

6.5.2 球面追踪(Sphere Tracing)

球面追踪是利用有符号距离场加速光线行进的技术。基本思想是：在每个位置，可以安全地前进到最近表面的距离。

算法流程：

pos = ray_origin
while (distance_traveled < max_distance):
    d = SDF(pos)
    if (d < epsilon):
        return hit
    pos += d * ray_direction
    distance_traveled += d
return miss

软阴影： 使用SDF可以高效计算软阴影： $$shadow = \min\left(1, k \cdot \min_{t \in [0,1]} \frac{SDF(p + t \cdot l)}{t}\right)$$ 其中 $k$ 控制阴影的软硬程度，$l$ 是到光源的方向。

环境光遮蔽(AO)： $$AO = 1 - \frac{k}{N} \sum_{i=1}^N \frac{SDF(p + s_i \cdot n)}{s_i}$$ 其中 $s_i$ 是采样距离，$n$ 是表面法向。

6.5.3 行军立方体(Marching Cubes)

行军立方体算法从体素数据（如Level Set场）提取三角网格。算法基于查表：立方体的8个顶点可以有256种内外配置。

基本步骤：

1. 分类：对每个体素，检查8个顶点的符号

cube_index = 0
for i in range(8):
    if values[i] < isolevel:
        cube_index |= (1 << i)

2. 边缘插值：在跨越等值面的边上线性插值 $$t = \frac{isolevel - v_1}{v_2 - v_1}$$ $$\mathbf{p} = \mathbf{p}_1 + t(\mathbf{p}_2 - \mathbf{p}_1)$$

3. 查表生成三角形：使用预计算的查找表

triangles = edge_table[cube_index]
vertices = []
for edge in triangles:
    vertices.append(interpolate(edge))

二义性处理： 某些配置存在二义性（如对角配置）。解决方法：

• 使用渐进立方体(Asymptotic Decider)
• 检查面和体的鞍点

改进算法：

• Dual Contouring：在对偶网格上生成顶点，保持尖锐特征
• Extended Marching Cubes：使用法向信息改善网格质量
• Adaptive Marching Cubes：根据曲率自适应细分

6.5.4 运动模糊与景深

这些效果模拟相机的物理特性，增强真实感。

运动模糊： 对快门开启时间积分： $$I(\mathbf{x}) = \frac{1}{T} \int_0^T L(\mathbf{x}(t), t) dt$$ 实现方法：

1. 后处理：使用速度缓冲区在屏幕空间模糊
2. 累积缓冲：渲染多个时间采样并平均
3. 随机采样：路径追踪时随机采样时间

时间采样的路径追踪：

t = random() * shutter_time
update_scene_to_time(t)
ray = generate_ray(pixel, t)
color += trace_ray(ray)

景深： 模拟有限孔径大小的透镜系统： $$I(\mathbf{x}) = \int_{aperture} L(\mathbf{x}, \mathbf{u}) K(\mathbf{u}) d\mathbf{u}$$ 其中 $K(\mathbf{u})$ 是孔径形状函数。

薄透镜模型：

lens_u, lens_v = sample_disk() * aperture_radius
ray_origin = camera_pos + (lens_u, lens_v, 0)
focal_point = ray_origin + focal_distance * initial_direction
ray_direction = normalize(focal_point - ray_origin)

散焦模糊的艺术控制：

• 光圈形状：圆形、多边形、自定义形状
• 散景效果：使用特殊的核函数
• 焦点呼吸：焦距随时间变化

6.6 体素渲染

体素渲染直接处理三维体数据，广泛应用于流体仿真、医学成像和科学可视化。与表面渲染不同，体渲染考虑整个体积内的贡献。

6.6.1 数字微分分析器(DDA)

DDA是三维空间中光线遍历体素的基础算法，类似于2D的Bresenham算法。目标是找出光线经过的所有体素。

3D DDA算法：

设光线方程为 $\mathbf{r}(t) = \mathbf{o} + t\mathbf{d}$，其中 $\mathbf{o}$ 是起点，$\mathbf{d}$ 是方向。

// 初始化
voxel = floor(ray_origin)
step = sign(ray_direction)
t_max = (voxel + (step > 0) - ray_origin) / ray_direction
t_delta = 1.0 / abs(ray_direction)

// 遍历
while (in_bounds(voxel)):
    visit(voxel)

    // 找到最小的t_max
    if (t_max.x < t_max.y && t_max.x < t_max.z):
        voxel.x += step.x
        t_max.x += t_delta.x
    elif (t_max.y < t_max.z):
        voxel.y += step.y
        t_max.y += t_delta.y
    else:
        voxel.z += step.z
        t_max.z += t_delta.z

优化技巧：

• 空间跳跃：使用层次结构跳过空体素
• 早期终止：当累积不透明度接近1时停止
• 自适应采样：根据内容密度调整步长

6.6.2 光线行进算法

光线行进是体渲染的核心，沿光线积分体素贡献。

基础算法：

color = 0
transmittance = 1.0
t = t_start

while (t < t_end && transmittance > epsilon):
    pos = ray_origin + t * ray_direction

    // 采样体素值
    density = sample_volume(pos)

    // 计算贡献
    dt = compute_step_size(t, density)
    absorption = exp(-density * dt)
    emission = compute_emission(density, pos)

    // 累积颜色
    color += transmittance * emission * (1 - absorption)
    transmittance *= absorption

    t += dt

自适应步长： 根据密度梯度调整步长： $$\Delta t = \min\left(\Delta t_{max}, \frac{\epsilon}{||\nabla\rho|| + \epsilon}\right)$$ 预积分传输函数： 预计算密度到颜色/不透明度的映射，存储为查找表： $$C(\rho) = \int_0^\rho c(s) \alpha(s) e^{-\int_0^s \alpha(u)du} ds$$

6.6.3 体积渲染方程

体积渲染考虑介质中的吸收、发射和散射。

辐射传输方程： $$\frac{dL}{ds} = -\sigma_t L + \sigma_a L_e + \sigma_s \int_{4\pi} p(\omega' \to \omega) L(\omega') d\omega'$$ 其中：

• $\sigma_t = \sigma_a + \sigma_s$：消光系数
• $\sigma_a$：吸收系数
• $\sigma_s$：散射系数
• $p(\omega' \to \omega)$：相函数
• $L_e$：发射

光学深度与透射率： 光学深度： $$\tau(s_0, s_1) = \int_{s_0}^{s_1} \sigma_t(s) ds$$ 透射率： $$T(s_0, s_1) = e^{-\tau(s_0, s_1)}$$ 体积渲染积分： 沿视线的最终辐射度： $$L = \int_0^D T(0, s) \cdot [\sigma_a(s)L_e(s) + \sigma_s(s)L_{in}(s)] ds$$ 单次散射近似： 忽略多次散射，只考虑直接光照： $$L_{in}(s) = L_{light} \cdot T(s, s_{light}) \cdot p(\omega_l \to \omega_v)$$

6.6.4 实时体渲染优化

实时应用需要各种优化技术来达到交互帧率。

空间跳跃(Empty Space Skipping)： 使用多分辨率数据结构标记空区域：

1. Min-Max块：存储每个块的最小/最大值

if (block_max < iso_min || block_min > iso_max):
    skip_block()

2. 距离图：预计算到最近非空体素的距离

skip_distance = distance_map[current_voxel]
t += skip_distance

GPU优化：

1. 纹理缓存：利用3D纹理硬件

value = texture3D(volume_texture, pos / volume_size)

2. 提前光线终止：

if (accumulated_opacity > 0.99):
    break

3. 块光线追踪：同时处理多条光线

for each 8x8 tile:
    find_entry_exit_points()
    march_rays_coherently()

时序一致性： 利用帧间相关性：

1. 重投影：使用上一帧结果作为初始猜测
2. 渐进式细化：低分辨率预览，逐步提高质量
3. 时间超采样：累积多帧结果降噪

混合渲染： 结合体渲染和表面渲染：

// 先渲染不透明表面
render_opaque_geometry()

// 然后渲染半透明体积
enable_blending()
for each volume:
    if (intersects_view_frustum(volume)):
        render_volume_with_depth_test()

这样可以正确处理体积与几何体的遮挡关系。

本章小结

本章深入探讨了高级输送格式和等势面方法，这些技术是构建高质量流体仿真系统的核心：

1. 高阶输送格式：WENO格式通过自适应权重选择避免数值振荡，TVD限制器保证数值稳定性，各种通量限制器在精度和稳定性间提供不同的权衡。

2. 有符号距离场(SDF)：提供了精确的几何表示，满足Eikonal方程 $||\nabla\phi|| = 1$，可以使用FMM或FSM高效构建，需要定期重新初始化以维持距离性质。

3. Level Set方法：将界面表示为零等势面，自然处理拓扑变化，但存在质量损失问题，可通过PLS等方法改善。

4. 自由表面追踪：VOF方法精确守恒但几何信息差，Level Set方法几何精度高但有质量损失，混合方法结合两者优势。

5. 渲染技术：路径追踪提供物理正确的渲染，球面追踪利用SDF加速，Marching Cubes从体数据提取网格，运动模糊和景深增强真实感。

6. 体素渲染：DDA算法高效遍历体素，光线行进沿路径积分，体积渲染方程考虑吸收、发射和散射，多种优化技术支持实时渲染。

关键公式回顾：

• Level Set演化方程：$\frac{\partial\phi}{\partial t} + V_n ||\nabla\phi|| = 0$
• 重新初始化方程：$\frac{\partial\phi}{\partial\tau} + S(\phi_0)(||\nabla\phi|| - 1) = 0$
• 渲染方程：$L_o = L_e + \int_{\Omega} f_r L_i (\omega_i \cdot \mathbf{n}) d\omega_i$
• 体积渲染积分：$L = \int_0^D T(0, s) \cdot [\sigma_a L_e + \sigma_s L_{in}] ds$

练习题

基础题

1. WENO权重计算 给定三个子模板的值：$u_{i-2} = 1$, $u_{i-1} = 2$, $u_i = 2.5$, $u_{i+1} = 2.8$, $u_{i+2} = 3$，计算五阶WENO格式在点$i+1/2$的数值通量。

提示：首先计算三个子模板的光滑度指标$\beta_k$，然后计算非线性权重$\omega_k$。

2. SDF性质证明 证明对于球面 $||\mathbf{x} - \mathbf{c}|| = R$，其有符号距离函数 $\phi(\mathbf{x}) = ||\mathbf{x} - \mathbf{c}|| - R$ 满足 $||\nabla\phi|| = 1$。

提示：直接计算梯度并求其模长。

3. Level Set质量计算 给定二维Level Set函数 $\phi(x,y) = x^2 + y^2 - 1$（单位圆），计算其所包围的面积。如果使用一阶精度的数值格式演化后，$\phi$ 变为 $\phi'(x,y) = x^2 + y^2 - 0.9$，计算质量损失百分比。

提示：面积为 $A = \int_{\phi<0} dxdy$。

4. 曲率计算 对于二维Level Set函数 $\phi(x,y) = x^2 + y^2 - R^2$，推导并计算界面上的曲率。

提示：使用公式 $\kappa = \nabla \cdot (\nabla\phi/||\nabla\phi||)$。

挑战题

5. 设计新的TVD限制器 设计一个新的通量限制器 $\phi(r)$，满足：

• TVD条件：$0 \leq \phi(r) \leq 2$
• 在光滑区域达到三阶精度
• 比van Leer限制器更少耗散

分析你的限制器的性质并与现有限制器比较。

提示：考虑分段多项式形式，确保在$r=1$附近的行为。

6. 优化FMM算法 快速行进法的标准实现使用堆来维护Trial集合。设计一种改进的数据结构，在网格规则且速度函数变化缓慢时能达到更好的性能。分析你的方法的时间复杂度。

提示：考虑利用网格的规则性和速度函数的局部性。

7. 实现PLS方法的粒子管理 设计一个高效的粒子管理系统用于PLS方法，需要处理：

• 粒子的空间索引（快速查找某区域内的粒子）
• 动态添加/删除粒子
• 并行化考虑

给出数据结构设计和关键算法的伪代码。

提示：考虑空间哈希或八叉树结构。

8. 体积渲染的重要性采样 对于参与介质的体积渲染，设计一种重要性采样策略，能够：

• 根据介质密度自适应采样
• 考虑光源位置进行方向采样
• 保证无偏性

推导你的采样概率密度函数和相应的权重。

提示：考虑使用光学深度进行距离采样，相函数进行方向采样。

常见陷阱与错误

数值格式陷阱

• WENO退化：当所有子模板的光滑度指标都很小时，权重可能退化，导致精度下降
• 限制器失效：在极值点附近，限制器可能过度限制，导致精度损失
• CFL条件违反：高阶格式可能需要更严格的CFL条件

SDF维护错误

• 重新初始化时机：过于频繁会累积误差，过于稀疏会失去距离性质
• 边界处理：FMM和FSM在边界附近需要特殊处理
• 数值精度：梯度计算的数值误差会累积

Level Set常见问题

• 质量损失累积：长时间演化后可能损失大量质量
• 薄结构消失：网格分辨率不足时，薄片、细丝等特征会消失
• 速度延拓错误：不当的延拓可能导致非物理的界面运动

渲染相关问题

• 采样不足：Monte Carlo方法的噪声，需要足够的样本数
• 数值精度：球面追踪的epsilon选择，过小导致无限循环，过大导致穿透
• 性能瓶颈：体渲染的内存带宽限制，需要优化内存访问模式

最佳实践检查清单

输送格式选择

• [ ] 评估问题的光滑性，选择合适阶数的格式
• [ ] 考虑计算成本与精度的平衡
• [ ] 验证格式在测试问题上的收敛阶
• [ ] 检查是否满足守恒性要求

SDF构建与维护

• [ ] 选择合适的初始化方法（解析/FMM/FSM）
• [ ] 设定合理的重新初始化频率和收敛准则
• [ ] 在界面附近使用更高的网格分辨率
• [ ] 验证梯度模长是否接近1

Level Set方法实施

• [ ] 使用高阶输送格式减少数值耗散
• [ ] 实施质量修正策略（如PLS）
• [ ] 正确处理速度场的延拓
• [ ] 监控质量守恒和界面锐利度

渲染优化

• [ ] 根据场景特点选择合适的渲染方法
• [ ] 实施空间加速结构
• [ ] 利用时间相关性和空间相关性
• [ ] 平衡质量与性能，使用自适应采样