摘要
本文提出一种流固耦合算法，能处理用三角面片表示的，无穷薄的固体。经典的流固耦合算法在三维格子上表示固体，但薄片没有内部区域，这就需要新方法。我们使用Robust ray casting来辅助插值，有限差分和渲染算法，使得液体不会穿过固体薄片。此外，我们提出了一种新方法，使得在enforcing incompressibility过程中，固体表面附近的流体不会被异常地压缩。这样，我们就可以模拟薄布料和少量水的互动。我们的方法既能处理刚体也能处理软体。我们提出了一种two-way coupling算法，让流体压力能影响到固体。 我们的算法既能用在规则网格上，也能用于自适应八叉树网格。
1 简介
略
2 先前研究综述
略
3 Robust One-Sided Interpolation of Data
我们需要阻止烟雾/水穿过薄软体。我们使用一种类似可见性检测的方法，用来确定哪些格点可用于插值和差分。可以使用robust ray casting算法做这一点，大体上就是把三角形面片增厚成楔形体，如[Bridson et al. 2002]所述，见图3左半边。取定插值点后，从该点做ray casting，把整个空间分成三个部分：可见部分、被遮挡部分、楔形体内部部分，如图3右半边所示。楔形体的可见点一定也是三角片的可见点，但在边界附近，三角片的可见点有可能是楔形体的内部点，甚至是被遮挡点。这相当于对物体边界做了模糊处理，因此能较为有效地消除舍入误差。在做插值和差分的时候，忽略所有被遮挡点和内部点，这就是单侧近似。
当固体移动的时候，某些点可能会被固体表面扫过，从一侧变到另一侧。在这种情况下，后续插值将把该点标为“无效”，因为这些点包含另一侧的信息。为了避免液体“漏过”表面，检测这些点至关重要，方法是每个三角形分别处理。在每一个时间步中，三角面片都会移动，在此过程中，如果一个格点穿过了三角形，就把它标记为“无效”。这种检测相当于求解一个三次方程，如[Bridson et al.2002]所述。为了健壮性，我们还需要把所有在时间步初始或结束时位于楔形体内部的点也标记为“无效”。所有时间步初始或结束时不在楔形体内，并且从三角片表面一侧变到另一侧的点，都被标记为和三角片发生碰撞。如果使用八叉树网格，分割出来的新节点，在插值中也被标为“无效”。而八叉树网格的合并操作不会新建节点，只会删除节点，所以没有特殊操作。
对于所有标为“无效”的点，我们都需要赋予它一个值，方法是使用Gauss-Jacobian迭代算法，把有效节点的值propagate出去。在每一次迭代中，每个“无效”节点都被赋予其周围一圈（一阶相邻）的可见有效节点的平均值，并被标记为“有效”。这种做平均的方法类似于其他人使用的混合方法，例如[Benson 1992].复杂的几何形状，或者折叠布料这样的操作，可能会让一些节点无论迭代几步也还是无效，因为它们没有任何相邻的可见有效节点。这时，我们需要再用Gauss-Jacobi迭代，不过在可见检测射线和固体相交时，需要特别地取值。例如，固体的速度，零密度，环境温度或固体温度，到固体的正距离等等。
我们使用标准的，和坐标轴对齐的，分层检测盒，用来加速相交检测。对于每个三角面片，我们都用一个稍大的bounding box框住它，用来把在流体模拟中，标记和它邻近的所有格子，这些格子可能需要特殊处理。
4 流体模拟
和[Fedkiw et al. 2001]提出的做法相同。忽略流体粘度。
虽然我们实现了八叉树网格上的算法，不过这里展示的算法主要用于规则网格。建议读者阅读[Losasso et al. 2004]获取细节。此外，我们实现了一个全新的，基于节点的流体求解器，包含速度、温度、密度、level set，所有值都定义在规则网格或八叉树网格的节点上。见4.2节。
4.1 计算中间速度场
我们使用半拉格朗日法，从[latex]\mathbf{x}[/latex]回溯到 [latex]\mathbf{x}-\Delta t\mathbf{u}[/latex]，用来计算节点速度。在这个过程中，使用两倍大小的楔形体进行单边插值（即图3中，[latex]\epsilon’=2\epsilon[/latex]）.这保证插值使用的点对于最终的速度节点一定是可见的。之后，向插值点周围8个格点发射的射线，就可以精确判断可见性（似乎是：结果点-插值点-周围格点，这样这样一个两段的折线）。如果某条射线和楔形体相交，就使用相交点处的固体速度。在中间速度场计算完成后，固体位置更新，把所有被固体移动过程中扫过的节点标记为“无效”，然后使用第3章中描述的方法，去propagate节点值给所有“无效”节点。
模拟水的时候，直接简单地给速度节点添加重力加速度。如果模拟烟，需要增加源项（因为烟会在温度带来的密度差作用下，向上飘散），取决于烟的密度[latex]\rho[/latex]和流体温度[latex]T[/latex]，也就是说，[latex]\mathbf{f}=-\alpha\rho\mathbf{z}+\beta(T-T_a)\mathbf{z}[/latex]，其中[latex]\mathbf{z}[/latex]指向上方，[latex]\alpha,\beta[/latex]是可调参数，[latex]T_a[/latex]是环境温度。烟的密度、流体的温度，都和速度场采用同样方式处理，如先前所述（对流步，可见性，相交检测，有效/无效节点，等等）。如果可见性检测射线和固体相交，就使用零密度和[latex]T_a[/latex]赋值。
为了在格点处计算vorticity confinement力（注：vorticity confinement是一种和投影法不同的流体模拟方法，不知道此处何意），我们需要使用六个相邻节点处的速度向量计算出涡度[latex]\omega=\nabla\times\mathbf{u}[/latex]。当然，如果某个节点不可见，就需要用相交点处的固体速度。然后就可以计算涡度的模长，再用六个相邻节点，做中心差分，计算模长的梯度[latex]\nabla|\omega|[/latex]。当然，和之前类似，如果节点不可见就用交点处固体的涡度。然后normalize这些梯度，得到vorticity location vector也就是[latex]\mathbf{N}[/latex]，从而在节点处计算涡度产生的源项[latex]\mathbf{f}=\hat{\epsilon}\Delta x(\mathbf{N}\times\omega)[/latex].
4.2 基于节点的流体求解器
我们在标准的MAC网格上求解压力。但是值都存在节点上，一种方法是在面中心做插值，然后求解压力，更新速度，最后把面上的速度场重新插值回节点。但是，额外插值会让求解质量变差，所以我们需要一种别的方法。我们在面上而非节点上初始化速度场。每一个时间步，我们都在计算中间速度场[latex]\mathbf{u}^*[/latex]之前把面上的速度场插值到节点上。然后，我们在对流步后并不会把中间速度场直接插值回去，而是计算速度差[latex]\Delta\mathbf{u}=\mathbf{u}^*-\mathbf{u}^n[/latex]，把这个修正增量插值回面上，然后进行更新（也就是FLIP）。这种方法比插值两次的方法要强得多。因为我们可以发现，来回插值会带来耗散，这种耗散和时间步无关，就算时间步长为零，耗散仍然会存在。而类似FLIP的增量法，计算出的增量则和时间步成比例，如果时间步趋于零就消失了。
这种节点/MAC混合求解法可以用于加速、简化任何NS方程求解器。不过在我们这里，还需要做几项改进。首先，在从面插值到节点时，需要对面做可见性检测，并在遮挡情况下使用固体上的速度。此外，在把更新速度场从节点插值回面上的时候，也需要对节点做类似的可见性检测。不过，如果有遮挡，就直接插值速度本身，而不是增量了，因为没法在固体上计算速度增量。
4.3 压力求解器
速度修正方程为
格子中心的压力值通过求解泊松方程得到：
离散化后这是一个对称矩阵，每个有流体的MAC格是其中的一行。如果流体是水，我们在水-空气的界面处设置Dirichlet边界条件[latex]p=1[/latex].而Neumann边界条件，意味着面中心处压强的导数为0，也就是说对速度不做更新。
如果没有正确处理边界条件，附着于固体上的薄层水膜，体积就很容易迅速变小，违反质量守恒。实际上，正确设置边界条件是保持质量守恒的最重要环节。因为如果边界上固体和流体速度不相等，那么流体就会流入或流出固体，导致质量不守恒。我们处理边界条件的方法，是论文的主要贡献之一。首先，我们发现，投影步得到的速度被用于下一个时间步，所以它应该和下一个时间步的固体速度相等。因此，我们先计算下一个流体时间步的步长，然后演化固体，固体时间步可能比流体的更短，所以可能有多个substep.演化完毕后，根据固体的位置更新，对每个固体节点计算出等效速度（位移除以流体时间步），然后把固体重新放回时间步开始时的位置。如此一来，这个等效速度就是固体在下一个时间步的速度，我们用它来设置Neumann边界条件，让边界处流体的速度等于这个速度。如此一来，薄层水膜的速度就能和固体速度完全一致，保证质量守恒。如图5所示，我们从格子中心向周围六个格子做射线，检测物体是否穿过两个压强节点的连线。如果穿过，我们就在这个面上设置Neumann边界条件，让流体速度等于之前计算的固体等效速度。然后按照标准流程计算散度，求解压力方程，并更新速度场。
在模拟很软的薄固体（例如布料在水流中运动）时，一个常见的困难，是固体会自己折叠起来，包住一部分流体，把它分离开来。可以很简单地以边界条件为起点，用flood fill算法检测出这种情况。如果一块流体被Neumann边界条件完全包裹，那么它求解压力用的矩阵，零空间将包含全1向量，从而不可逆。有一种共轭梯度法可以用来求解此类矩阵，不过在此之前，需要先施加一个兼容性条件[Peyret and Taylor 1983].这个操作需要对每个有零空间的区域都做一次，使用边界格子面上的速度，计算出该区域的单位面积净流入和净流出。然后对于每个边界格子面，使用这个值和这个面的面积，计算出新的临时速度，保证穿过区域边界的净流量为零。最后，求解压力方程，保证这个区域散度为零。
4.4 水
我们使用particle level set [Enright et al.2002]方法模拟水。[latex]\phi\leq 0[/latex]是水，[latex]\phi>0[/latex]是空气。对于水，只需要求解速度场。每个时间步，我们都会把节点上的速度向表面之外，外推3~5个格子，以得到速度边界条件。做法是：先把这个3~5格宽的窄带中，所有格子按照[latex]\phi[/latex]值从小到大排序（当然，这需要在level set的reinitialization之后做）。然后按此顺序，逐个节点求解方程[latex]\nabla\phi\cdot\nabla\mathbf{u}=0[/latex].为了避免流体渗漏，我们排除掉所有不可见的边界值。有些节点可能没有可见的邻居，此时把它们也标为不可见。在外推完毕之后，再给所有不可见节点赋值，方法见上面第3章。
4.4.1 Level Set Method
根据[Enright et al.2005]，particle level set方法靠粒子提供精确的距离函数值，靠level set提供连通、平滑的表面。他们在文中称，可以把level set的高阶advection换成semi Lagrangian，而不影响精度。Level set的值存在节点上，因此semi Lagrangian和速度场advection类似。在用周围8个节点值（也就是正方体的8个顶点）进行插值的时候，如果碰到不可见节点，就用其他7个点的值计算出不可见节点的值。方法是，先看8个顶点中，那个不可见节点的一阶相邻点（有3个），是否对插值点可见。如果可见，就取它们的平均值，作为不可见节点的值。否则，就检查二阶相邻点（有3个），再不行就检查三阶相邻的（有1个）。如果都不行，就意味着正方体里面没有任何一个点对插值点可见，那么就不应该更新原先那个节点的值。在这种情况下，我们把那个节点的[latex]\phi[/latex]置为无效，再单独处理（见第3章）。
Level set储存的值是有向距离函数，使用fast marching进行维护（[Osher and Fedkis 2002]）.通过检测相邻格子的符号变化，就可以知道哪些格子在表面附近。如果一个格子[latex]\phi\leq 0[/latex]，且有一个不可见的相邻格，我们也认为它在表面附近。但如果[latex]\phi>0[/latex]，且没有一个可见相邻格的 [latex]\phi\leq 0[/latex] ，就不算是表面附近。这两条特判，让我们得以考虑流体-固体的表面，并且保证流体不会穿透固体，影响外面的空气。我们给所有表面附近的格子赋[latex]\phi[/latex]的初值，依据是它在和相邻格子符号相反的那个坐标轴上，离表面有多远。如果当前节点[latex]\phi>0[/latex]，我们应当忽略掉相邻格子不可见的坐标轴。如果当前节点[latex]\phi\leq 0[/latex]，我们取“离固体的距离”和“不可见方向上的[latex]|\phi|[/latex]”的最小值。这个特判，可以避免水“粘在”固体上。在level set初始化后，我们用fast marching算法，把它推广到整个计算域。不过在更新某点的level set值时，需要排除掉对它不可见的相邻节点（和速度场的extrapolation类似）。
4.4.2 Particle Level Set Method
[latex]\phi[/latex]值为负的粒子（也就是水）需要和固体进行碰撞检测，用来避免水漏过固体。这里我们使用碰撞距离法：给每个粒子预先设定一个碰撞距离值，是一个[latex]0.1\Delta x[/latex]到 [latex]\Delta x[/latex] 之间的随机数。在对某个粒子做碰撞检测的时候，我们先在固体上找到离这个粒子最近的点，并计算出固体在该点的法向。如果粒子距离固体的距离，小于这个粒子的碰撞距离，那么我们就把它沿着法向，向外移动，直到达到碰撞距离。如果在移动过程中，粒子和某个固体相交，就删除这个粒子，或不移动它。我们可以发现，正确实现“负粒子”和固体的碰撞检测，可以大大增加求解固体表面薄层水膜的能力。
粒子的速度通过插值决定，需要像4.1节描述的那样，对相应的8个节点发出射线，做可见性检测。对于离固体距离小于碰撞距离的“负粒子”，我们应当直接调整其速度的法向分量，使之不会继续靠近固体。[Enright et al. 2005]指出，二阶精度的粒子演化算法是很重要的，在用semi-Lagrangian做advection的情况下尤其如此。因此，我们先正向演化粒子，期间和静止的固体做碰撞检测。然后，我们用原先的速度场计算粒子在新位置的速度，再把粒子移回原先位置，这样精度就是二阶的。对于“负粒子”，如果它在新位置或老位置距离固体太近，低于碰撞距离，就应该调整其速度的法向分量。
在执行如上所述的，二阶精度粒子演化步之前，我们需要检查粒子和固体移动轨迹是否相交，删除掉所有相交的正粒子，但尝试根据相交的那个三角面片，去调整负粒子。根据粒子和面片的初始位置，我们记下粒子在三角的哪一侧。然后根据粒子和面片的终止位置，我们沿三角面片的法向去移动粒子，让它呆在同一侧，并且沿法向离三角面片的距离，应当超过碰撞距离。最后，我们检查粒子新的移动路径，看它是否和移动的固体相交，如果相交，就删去这个粒子。在advection过后，应当调整所有的负粒子位置，使它们距离固体至少有碰撞距离那么远，方法如前所述。
在更新level set和粒子之后，我们使用粒子上携带的值更新level set值。这里需要考虑碰撞，并且只使用可见的粒子。最后，应当根据[latex]\phi[/latex]的值调整粒子的半径，并且删去离表面太远的粒子。方法也是在粒子中心处计算[latex]\phi[/latex]值，和semi-Lagrangian射线类似。每模拟10~20帧，就需要reseed粒子，用来更好地表示表面。首先不管固体，直接reseed，然后计算每个粒子中心处的[latex]\phi[/latex]，删去符号错误的粒子。
5 布料和薄壳模拟
本文不提出特殊的固体模拟方法，而是把它视作“黑箱”。因此，你可以任意采用自己喜欢的固体模拟技术，这和流体模拟、耦合部分独立。固体模拟算法只需要给出在若干离散时间点（也就是每一个时间步）处，固体mesh上每个节点的位置，由此就可以计算出每个固体节点的速度。如果需要计算三角面片内部某点的速度，就使用重心坐标系进行插值。对于刚体模拟，我们使用[Guendelman et al. 2003]的方法。我们在文中展示的例子，固体模拟算法不超过[Hahn 1988; Moore and Wilhelms 1988].对于布料模拟，我们使用[Bridson et al. 2003]提出的布料模型，包括其屈服公式，和[Bridson et al. 2002]提出的自碰撞算法。
5.1 流固耦合
传统的流固耦合算法，是把固体速度作为流体速度的边界条件，然后把流体速度作为固体的边界条件[Benson 1992].正如4.3节中所述，确保流体和固体移动速度相同很重要，只有这样才能正确模拟薄层水膜，不致损失体积。对于可压缩流体，常见方法是在固体表面上施加流体的压力，但对于不可压缩流体，[Fedkiw 2002]指出，这个压力太僵硬，噪声太大。因此我们提出一个新算法，用来获取薄片状刚体或软体受到的流体力。
如4.3节所述，我们从每个格子中心向相邻六个格子中心做射线，用来检测固体是否穿过两个压强相邻节点的连线，如图5所示。如果穿过，我们就把相交点的固体速度值赋给两个格子中间的面。此外，我们还用固体密度乘以这个面的面积，计算出这个面的质量m。然后，我们用m除以这个面的控制体积，也就是它两个相邻格子体积之和的一半（如果用八叉树，需要除以恰当的系数）。对于规则网格，就是[latex]V=\Delta x\Delta y \Delta z[/latex]了。最后，把这个面的密度设为[latex]m/V[/latex].其他一般的面，密度就是流体密度。由于流体的密度可变，锁门需要求解一个变系数泊松方程，如[Nguyen et al. 2002]所述，其中相应方程为：
我们不像[Carlson et al. 2004]和[Zhu and Peskin 2002]那样，用这里的流体速度去移动固体。这样的主要优点，是我们不需要研究如何计算刚体对浮力、碰撞、弹性之类的反应，而是只需要简单地用“黑箱”模拟固体。此外，我们的固体不受流体网格分辨率的限制，网格粗糙一点，只是导致流体向固体施加的力比较“糊”，但仍然可以把固体本身算得很精细。
在求解公式6之后 ，为了计算固体受力，我们首先对每个是固体一部分的面，计算出其上的压差。然后我们把压差从面上转移到节点上，方法是在每个节点处，把与之相邻的至多四个面上的值做平均。之后，我们使用Gauss-Jacobi平均化算法，类似于之前给无效节点赋值那样，把这个压差扩散到固体附近窄带上的格点。这个过程对于三个维度各自独立做一遍。然后对于每个三角面片，在其重心处，线性插值出压差，把压差矢量和三角形法向做点积，就得到了这个面片上的压强。最后，用这个压强乘以面片面积，再乘以面片法向，就得到了面片受到的压力值。对于刚体，就可以把所有面片累加起来，得到总受力和总力矩。对于软体，就需要把面片受到的力均匀分配给三个节点。
在求解公式4的时候，需要在整个薄片状固体上施加Neumann边界条件，所以计算出来的压力可能会有显著跳变和噪声，原因是压力必须迫使流体的速度和固体一致，保持质量守恒。但在求解公式6的时候，并不施加Neumann边界条件，压强值按照公式(5)去修正流体速度，使之和固体一致（虽然我们不直接这么修正固体的速度）。在求解椭圆方程形式的压力方程时，我们能够考虑到流体对固体的影响，这让我们的耦合算法更加稳定，并且让压力场更加平滑。
6 算法总结
整个算法如下。
第一步，使用固体在时间步n和n+1处的位置信息，计算处中间速度场[latex]\mathbf{u}^*[/latex].然后把 [latex]\mathbf{u}^*[/latex] 代入公式6中的[latex]\mathbf{u}^{old}[/latex]，耦合求解压强场。使用这个压强场，计算出固体的受力，据此把固体从第n+1步演化到第n+2步，根据n+1，n+2步固体的位置，计算出它的等效速度。使用这个等效速度，对流体的中间速度场做投影，使之满足无散条件。
然后重复上面这个过程。