GAMES103课程笔记10-Surface Waves

这个系列是GAMES103：基于物理的计算机动画入门(GAMES 103: Intro to Physics Based Animation)的同步课程笔记。本课程是计算机动画的入门课程，着重介绍各种基于物理的动画仿真模拟技术。本节课介绍物理仿真中的水波模拟。

Fluid Effects

流体仿真(fluid simulation)是物理仿真中重要的研究内容。现实生活中常见的流体包括液体、气体、烟尘等。和前面介绍过的固体相比，流体具有更加丰富的行为，因此流体仿真也往往要更加复杂。

根据认识问题的角度不同，流体仿真的方法大体可以分为Lagrangian方法(Lagrangian approach)以及Eulerian方法(Eulerian approach)两大类。简单来说，Lagrangian方法更多地是从粒子的角度来模拟流体的行为，即考虑每个粒子在流场中的运动；而Eulerian方法则更多地从空间的角度来模拟流体运动，对于固定的空间Eulerian方法会考虑流场的不同物理量在空间中的变化。尽管本节课不会涉及这两类方法的细节，但在后面的课程中则会深入讨论二者的差异。

Height Field Model

高度场(height field)是进行水波模拟的经典方法。以2D(1.5D)空间为例，我们可以把水波看作是沿 $x$ 轴分布的高度函数 $h (x)$ 。类似地，水波的速度也是一个一维函数 $u (x)$ 。

显然只要我们计算出每一时刻的高度场和速度场就可以实现对水波的模拟。对于高度场，在任意位置和任意时刻需要满足微分方程：

\frac{d h (x)}{d t} + \frac{d (h (x) u (x))}{d x} = 0

为了便于理解，我们可以对上式进行变形：

d h (x) d x = - d (h (x) u (x)) d t

等号左边表示任意时刻 $x$ 位置流量的变化量。而右边 $d (h (x) u (x))$ 项则可以表示为：

d (h (x) u (x)) = h (x + d x) u (x + d x) - h (x) u (x)

即 $x$ 位置单位时间的流出量，对它取负号并乘上 $d t$ 即为净变化量，因此左右两端相等。

而对于速度场则需要满足方程：

\frac{d u (x)}{d t} = - u (x) \frac{d u (x)}{d x} - \frac{1}{ρ} \frac{d P (x)}{d x} + a (x)

实际上上式即为无黏性流体的一维NS方程，其中 $a (x)$ 表示外力。

Shallow Wave Equation

对于高度场方程，我们将 $d (h u)$ 展开可以得到：

\frac{d h}{d t} + \frac{d (h u)}{d x} = \frac{d h}{d t} + h \frac{d u}{d x} + u \frac{d h}{d x} = 0

对于浅水波(shallow wave)我们认为在高度 $h$ 在局部的变化量为0，得到化简后的方程：

\frac{d h}{d t} + h \frac{d u}{d x} = 0

对于速度场方程，我们忽略advection和外力可以得到化简形式：

\frac{d u}{d t} = - \frac{1}{ρ} \frac{d P}{d x}

将两个方程分别对 $t$ 和 $x$ 求导，得到：

\frac{d^{2} h}{d t^{2}} + h \frac{d^{2} u}{d x d t} = 0

\frac{d^{2} u}{d x d t} = - \frac{1}{ρ} \frac{d^{2} P}{d x^{2}}

消去公共项 $\frac{d^{2} u}{d x d t}$ 就得到了shallow wave equation(SWE)：

\frac{d^{2} h}{d t^{2}} = \frac{h}{ρ} \frac{d^{2} P}{d x^{2}}

Discretization

通过求解SWE就可以实现对水波的仿真。在计算导数时我们需要对自变量进行离散，具体地对于任意一元函数 $f (t)$ 我们可以通过差分的方式来计算导数的近似值。

对于二阶导项，我们则需要再对一阶导进行差分：

这样经过整理最终得到离散化的SWE：

h_{i} (t_{0} + Δ t) = 2 h_{i} (t_{0}) - h_{i} (t_{0} - Δ t) + \frac{Δ t^{2} h_{i}}{Δ x^{2} ρ} (P_{i + 1} + P_{i - 1} - 2 P_{i})

Volume Preservation

除了SWE外在仿真时还需要注意体积守恒(volume preservation)，即任意时刻流体的总体积不能增加或减少。我们假定 $t$ 时刻和 $t - Δ t$ 时刻流体的体积为 $\sum_{i} h_{i} (t) = \sum_{i} h_{i} (t - Δ t) = V$ ，根据SWE可以计算出 $t + Δ t$ 时刻的体积：

\begin{aligned} \sum_{i} h_{i} (t + Δ t) & = 2 \sum_{i} h_{i} (t_{0}) - \sum_{i} h_{i} (t_{0} - Δ t) + \sum_{i} \frac{Δ t^{2} h_{i}}{Δ x^{2} ρ} (P_{i + 1} + P_{i - 1} - 2 P_{i}) \\ = V + \sum_{i} \frac{Δ t^{2} h_{i}}{Δ x^{2} ρ} (P_{i + 1} + P_{i - 1} - 2 P_{i}) \end{aligned}

上式说明仿真时流体的体积会随第二项 $\sum_{i} \frac{Δ t^{2} h_{i}}{Δ x^{2} ρ} (P_{i + 1} + P_{i - 1} - 2 P_{i})$ 在 $V$ 附近波动，要满足体积守恒则需要保证第二项恒为0。

要保证体积守恒的一种办法是保证相邻格子上流体的流入和流出量相等，因此需要稍微修改SWE：

h_{i} (t_{0} + Δ t) = 2 h_{i} (t_{0}) - h_{i} (t_{0} - Δ t) + \frac{Δ t^{2}}{Δ x^{2} ρ} (\frac{h_{i - 1} + h_{i}}{2} (P_{i - 1} - P_{i}) + \frac{h_{i + 1} + h_{i}}{2} (P_{i + 1} - P_{i}))

实践中更常用的方法是使用一个固定的常量 $H$ 来代替 $h_{i}$ ，这样也能够保证更新后流体总体积不变：

h_{i} (t_{0} + Δ t) = 2 h_{i} (t_{0}) - h_{i} (t_{0} - Δ t) + \frac{Δ t^{2} H}{Δ x^{2} ρ} (P_{i + 1} + P_{i - 1} - 2 P_{i})

Pressure

对于压力项 $P_{i}$ ，我们可以利用压强公式 $P = ρ g h$ 来计算：

\begin{aligned} h_{i} (t_{0} + Δ t) & = 2 h_{i} (t_{0}) - h_{i} (t_{0} - Δ t) + \frac{Δ t^{2} H}{Δ x^{2} ρ} (P_{i + 1} + P_{i - 1} - 2 P_{i}) \\ = 2 h_{i} (t_{0}) - h_{i} (t_{0} - Δ t) + \frac{Δ t^{2} H g}{Δ x^{2}} (h_{i + 1} (t_{0}) + h_{i - 1} (t_{0}) - 2 h_{i} (t_{0})) \\ = 2 h_{i} (t_{0}) - h_{i} (t_{0} - Δ t) + α (h_{i + 1} (t_{0}) + h_{i - 1} (t_{0}) - 2 h_{i} (t_{0})) \end{aligned}

Viscosity

如果要模拟流体的黏性则可以在更新时为动量项添加阻尼：

\begin{aligned} h_{i} (t_{0} + Δ t) & = h_{i} (t_{0}) + (h_{i} (t_{0}) - h_{i} (t_{0} - Δ t)) + α (h_{i + 1} (t_{0}) + h_{i - 1} (t_{0}) - 2 h_{i} (t_{0})) \\ \to h_{i} (t_{0}) + β (h_{i} (t_{0}) - h_{i} (t_{0} - Δ t)) + α (h_{i + 1} (t_{0}) + h_{i - 1} (t_{0}) - 2 h_{i} (t_{0})) \end{aligned}

将前面的步骤综合起来就可以得到水波仿真算法：

Boundary Conditions

要求解SWE还需要考虑边界条件。一般来说边界条件有两种类型：Dirichlet边界和Neumann边界。Dirichlet边界是指已知函数在边界处的值，而Neumann边界则是已知函数在边界处的微分。对于Dirichlet边界条件，它等价于在边界处固定了外部流体的液面高度，因此适用于模拟开放水面；而对于Neumann边界一般情况下我们会假定边界处的微分为0，即没有其他水体流入，因此Neumann边界适用于模拟封闭水面。