这个系列是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 + dx) u(x + dx) - h(x) u(x)\]即\(x\)位置单位时间的流出量，对它取负号并乘上\(dt\)即为净变化量，因此左右两端相等。
而对于速度场则需要满足方程：
\[\frac{d u(x)}{d t} = -u(x) \frac{d u(x)}{d x} - \frac{1}{\rho} \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}{\rho} \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}{\rho} \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}{\rho} \frac{d^2 P}{d x^2}\]Discretization
通过求解SWE就可以实现对水波的仿真。在计算导数时我们需要对自变量进行离散，具体地对于任意一元函数\(f(t)\)我们可以通过差分的方式来计算导数的近似值。
对于二阶导项，我们则需要再对一阶导进行差分：
这样经过整理最终得到离散化的SWE：
\[h_i(t_0 + \Delta t) = 2 h_i(t_0) - h_i (t_0 - \Delta t) + \frac{\Delta t^2 h_i}{\Delta x^2 \rho} (P_{i+1} + P_{i-1} - 2 P_i)\]Volume Preservation
除了SWE外在仿真时还需要注意体积守恒(volume preservation)，即任意时刻流体的总体积不能增加或减少。我们假定\(t\)时刻和\(t - \Delta t\)时刻流体的体积为\(\sum_i h_i(t) = \sum_i h_i(t - \Delta t) = V\)，根据SWE可以计算出\(t + \Delta t\)时刻的体积：
\[\begin{aligned} \sum_i h_i(t + \Delta t) &= 2 \sum_i h_i(t_0) - \sum_i h_i (t_0 - \Delta t) + \sum_i \frac{\Delta t^2 h_i}{\Delta x^2 \rho} (P_{i+1} + P_{i-1} - 2 P_i) \\ &= V + \sum_i \frac{\Delta t^2 h_i}{\Delta x^2 \rho} (P_{i+1} + P_{i-1} - 2 P_i) \end{aligned}\]上式说明仿真时流体的体积会随第二项\(\sum_i \frac{\Delta t^2 h_i}{\Delta x^2 \rho} (P_{i+1} + P_{i-1} - 2 P_i)\)在\(V\)附近波动，要满足体积守恒则需要保证第二项恒为0。
要保证体积守恒的一种办法是保证相邻格子上流体的流入和流出量相等，因此需要稍微修改SWE：
\[h_i(t_0 + \Delta t) = 2 h_i(t_0) - h_i (t_0 - \Delta t) + \frac{\Delta t^2}{\Delta x^2 \rho} \bigg( \frac{h_{i-1} + h_i}{2} (P_{i-1} - P_i) + \frac{h_{i+1} + h_i}{2} (P_{i+1} - P_i) \bigg)\]
实践中更常用的方法是使用一个固定的常量\(H\)来代替\(h_i\)，这样也能够保证更新后流体总体积不变：
\[h_i(t_0 + \Delta t) = 2 h_i(t_0) - h_i (t_0 - \Delta t) + \frac{\Delta t^2 H}{\Delta x^2 \rho} (P_{i+1} + P_{i-1} - 2 P_i)\]
Pressure
对于压力项\(P_i\)，我们可以利用压强公式\(P = \rho g h\)来计算：
\[\begin{aligned} h_i(t_0 + \Delta t) &= 2 h_i(t_0) - h_i (t_0 - \Delta t) + \frac{\Delta t^2 H}{\Delta x^2 \rho} (P_{i+1} + P_{i-1} - 2 P_i) \\ &= 2 h_i(t_0) - h_i (t_0 - \Delta t) + \frac{\Delta t^2 H g}{\Delta 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 - \Delta t) + \alpha (h_{i+1}(t_0) + h_{i-1}(t_0) - 2 h_i(t_0)) \end{aligned}\]Viscosity
如果要模拟流体的黏性则可以在更新时为动量项添加阻尼：
\[\begin{aligned} h_i(t_0 + \Delta t) &= h_i(t_0) + (h_i(t_0) - h_i (t_0 - \Delta t)) + \alpha (h_{i+1}(t_0) + h_{i-1}(t_0) - 2 h_i(t_0)) \\ &\rightarrow h_i(t_0) + \beta (h_i(t_0) - h_i (t_0 - \Delta t)) + \alpha (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边界适用于模拟封闭水面。
对于Neumann边界为0的情况我们只需要跳过边界即可，相应的仿真算法如下：
Coupling
当水面上存在其它物体时会产生耦合(coupling)的现象，即物体会排开水体产生波浪，而波浪又会影响物体的运动状态。
对于流体而言，我们需要计算出物体排出的体积并把这部分流体排出去。要模拟这种排开水的行为我们可以假定在物体存在的位置有一个虚拟的液面高度\(v\)，假设初始液面高度为\(h\)并记排开液体的高度为\(e\)，根据SWE可以得到方程组：
整理方程可以得到线性方程组(Poisson方程)：
通过求解Poisson方程即可实现流体部分的仿真，对应的流程如下：
而对于固体部分，只需要把排开液体的重量作为浮力施加到相应位置即可：
