作者： wmdcstdio

传给glutMouseFunc()的(x,y)是左上至右下的坐标：x表示离屏幕左边缘的像素数，y表示离屏幕上边缘的像素数，如图：

而glReadPixels()采用的坐标系统则与此不同，它是从左下到右上的。也就是说，x表示离屏幕左边缘的像素数，y表示离屏幕下边缘的像素数，如图：

所以，如果想把glutMouseFunc()取得的坐标传给glReadPixels()，就需要做一个变换：

y=Height-y-1

3 SOAP FILM DYNAMICS WITH VARYING THICKNESS

3.1 Definitions and Setting

流体在薄膜上的流动遵循带有表面张力项的无粘NS方程：

这个x上面一点就是速度的意思。(1)右端第一项是表面张力。里面的u代表气-液界面，l代表液-气界面。H是平均曲率，n是法向。这里u和l得名于“upper”和“lower”.

然后[latex]\delta_{S^u},\delta_{S^l}[/latex]是狄拉克函数，在[latex]S^u,S^l[/latex]之外消逝。

压强p在除边界之外的地方二阶可导，在每个泡泡里面的空气区域中取恒定值。 S是surface：

这里面A是[latex]R^2[/latex]的子集。然后在一个点可以确定局部坐标系：

这里[latex]S_a[/latex]和[latex]S_b[/latex]代表偏导。

厚度场为h，因此薄膜上下表面为：

然后它定义了一个[latex]\Gamma[/latex]：

这个[latex]\Gamma[/latex]实际上是一个“柱”，也就是以S(a,b)表示的这个表面点为中心，向上下到薄膜表面的一个线段。

然后X就是这个“柱”的重心：

文中将用X代表表面S。

3.2 Normal Acceleration

因为厚度h足够小，所以可以忽略其上的变化，用表面上点的速度代表这个“柱”的速度：

单看法向加速度的话，根据NS方程，就是：

最后一个等号是根据狄拉克函数的性质，以及如下等式：

这里面，[latex]\nabla^2_s[/latex]是表面S上的Laplace-Beltrami算子。我们认为薄膜两边的body force一样，因为它足够薄。然后也假设[latex]\nabla_{S^u}^2S^u+\nabla_{S^l}^2S^l=2\nabla_S^2S[/latex].

然后考虑压强项。在薄膜表面，也就是“柱”[latex]\Gamma[/latex]的两端，根据Ishida et al’s [2017]的近似：

也就是说，压强差乘以法向。然后把中间的压强z方向梯度积分起来，其实也等于两端压强差：

我们认为法向在整个“柱”上都是一样的，也就是[latex]n^u=n^l=n[/latex]，结果就是：

其实就是两端压强差，这么点事……

这样的话，合起来的加速度公式就是：

影响这个法向运动的主要是：表面张力（曲率）、两端压强差、外力。

3.3 Tangential Acceleration

切向加速度就是：

这里面我们认为压强p是二阶可导的，同样[latex]n^u=n^l=n[/latex].

然后由于法向的x，y分量为零，也就是：

这里面的[latex]\nabla_S[/latex]是表面梯度（surface gradient），P是一个“柱”的压强积分：

注意：切向没有表面张力。也就是说，文章只算毛细力，忽略了Marangoni力。

3.4 Thickness Evolution

首先回顾厚度的定义：

然后我们看，在[latex]x^u[/latex]处的法向力为：

类似地

把这两个合到一起算出来h的变化率：

这里把上下表面的表面梯度之差写成了[latex]\nabla_S^2h[/latex]。然后

可以发现，如果忽略掉Q，h的方程就是一个advected wave equation.然后Q的作用是“挤压”这个film，保持不可压性，就像流体求解器里面的压强场一样。

3.5 Pressure

压差[latex]\nabla p[/latex]实际上是每个泡泡体积守恒产生的一个修正项。我们使用Ishida et al. [2017]的方法做这个。

然后我们看方程(20)里的P。根据质量守恒，厚度h满足：

使用一个类似投影法的方法。首先对(20)做对流，然后加上外力，然后剩下：

这里带波浪线的是对流&加外力之后的速度场。把它和(26)联立可以得到：

这个就相当于N-S求解器里的投影了。这个和拿triangle mesh做的二维不可压流体求解器几乎一致，但是它允许质量在切向速度和法向的h变化率之间转意。所以，速度场本身是可压的，散度会被h的变化率吸收掉。

本文的算法没有真正地计算Q。它是用Q去修正数值误差导致的体积变化。每个时间步结束之后，用Q去scale up厚度项h，保证薄膜的质量守恒。

3.6 Evaporation and Bursting Bubbles

蒸发就是每步把h削薄一点：

破裂就是当某个地方的h太小的时候，删掉一块表面。把这个bubble的质量分给周围其他bubble.

4 IMPLEMENTATION

大致可以总结为：

1、切向advection.
2、根据切向速度和h算P。
3、根据算出来的P更新切向速度。
4、Advect厚度。
5、给h加表面张力。
6、用Q对h做体积修正。
7、计算法向速度，其中有表面张力，需要用到h。
8、做空气体积守恒。

7 约束和碰撞

我们使用点约束（point constraints），这样一方面施加软体约束（例如bone attachments），一方面解决物体之间的碰撞/自碰撞。特别地，我们在规则网格（lattices）上放一些代理点（proxy points），在这些点上记录力，并用三次插值法把这些力分配回六面体网格节点上（hexahedral cells）。[Sifakis et al. 2007]展示了这种方法的效率。

碰撞检测

对碰撞的响应（collision response）由一系列代理点决定，它们大致覆盖了碰撞表面。我们使用一种基于惩罚项（penalty）的响应方式，这个惩罚项取决于代理点处的穿透深度，以及该点的法向。对于刚体而言，我们直接用level set表示它，就能计算出代理点处的穿透深度。对于自碰撞（self-collision），由于物体迅速形变，导致难以每一步都精确地算出一个有向距离函数（signed distance function）.

自碰撞穿透深度计算

在每个代理点处，我们先计算出来，它目前所处的那个，变形后的六面体格子。这一步可以用一个沿坐标轴对齐的（axis aligned）包围盒（bounding box）树结构来完成，其每个叶子节点包含一个变形后的六面体格。为了避免假阳性，我们不在物体坐标系（material coordinates）下去找周围一圈27个六面体格。在代理点附近的每个六面体格，都被切分成四面体（tetrahedralized），以计算出重心坐标，确定代理点在物体上的位置。对于每个物质点，我们查询之前储存的，静息状态下的level set值[latex]\phi_0[/latex].如果有多个负值，我们就取其中最接近零的那个[latex]\phi_0[/latex]，来计算最接近静息状态下，表面上与之最接近的点。然后我们去找这个表面上的点，在变形后的位置，用来计算穿透深度[latex]|x_s-x_p|[/latex]，其中[latex]x_s[/latex]是变形后的表面最近点，[latex]x_p[/latex]是代理点。

碰撞响应

对于自碰撞和互碰撞两种场景，我们都采用同样的方式处理。即：“生成”一个静息长度（rest length）为零的弹簧，一头连接代理点，一头连接表面最近点。其杨氏模量可以是各向异性的，突出法向。特别地，其能量函数为：[latex]\Psi(x_p,x_s)=k(x_p-x_s)^TM(x_p-x_s)/2[/latex]，其中[latex]M=(1-\alpha)nn^T+\alpha I[/latex]，其中[latex]\alpha\in[0,1][/latex].如果[latex]\alpha=1[/latex]那就是传统的各向同性弹簧。如果[latex]\alpha=0[/latex]那么就是标准的，点斥力（point repulsion）.这种stiffness的各向异性，允许物体沿平面运动（垂直于穿透方向）。在实践中，我们发现[latex]\alpha\in[0.1,0.5][/latex]对自碰撞效果最好，而[latex]\alpha=0[/latex]对物体之间的互碰撞就足够了。

1 模拟基础

用“状态向量”表示一个质点的状态：

如果三维，就是六个数。这个Y可以扩展到n个质点：

我们先暂且考虑一个质点。F(t)是该质点t时刻所受合力。假设该点有质量m，那么状态向量的变化率就是：

2 刚体基本概念

2.1 位置和方向

质点没有方向，刚体有，所以复杂一点。描述刚体，首先需要位置x(t)，然后是3*3旋转矩阵R(t).我们称x(t)和R(t)是刚体的“空间变量”。

在固定的body space里面定义刚体的形状。之后就可以用x(t)和R(t)把它变换到world space下。我们约定刚体质心在原点(0,0,0).假设R(t)描述了刚体关于质心的转动。那么，body space下的位置r，将被变换到R(t)r的位置。如果p0是body space里任意一点，那么它在world space的位置就是

R的物理意义是什么呢？考虑body space下的x轴，也就是(1,0,0).它在t时刻world space下的朝向是：

把R的分量写作

那就是

当然y和z轴被变换到的位置就是

也就是说，R矩阵的三列，分别是world space下，初始时body space的三根坐标轴。

2.2 线速度

2.3 角速度

用矢量[latex]\omega(t)[/latex]表示角速度。其朝向为转动轴，长度表示转动速率。把刚体上的点r分解成和转动轴[latex]\omega(t)[/latex]平行的a，以及与之垂直的b。假设刚体不平动，那么r就在绕轴转动，这个转动的线速率是[latex]|b||\omega(t)|[/latex]。由于b和[latex]\omega[/latex]垂直，所以有

也就是说[latex]\dot{r}(t)=\omega(t)\times(b)[/latex].然后[latex]\omega[/latex]叉乘a得零（因为平行），所以：

所以转动导致的线速度就是

根据这个我们发现，R第一列（也就是x轴）的变化率是：

其他两个轴也一样，这就是

这个写法太麻烦，下面简化一下。根据叉乘公式，[latex]a\times b[/latex]结果为

对于矢量a，定义矩阵a*：

那么

那么[latex]\dot{R}(t)[/latex]就是

即

也就是

这个形式写起来比较优雅：[latex]\dot{r}(t)=\omega(t)\times r(t)[/latex]，[latex]\dot{R}(t)=\omega(t)*R(t)[/latex]

2.4 刚体质量

假设刚体由一堆质点构成。i号质点在world space下的位置是

当然这个[latex]r_{0i}[/latex]是它在body space的位置。

然后总质量就是

2.5 质点的速度

对(2-13)求时间导数，注意到刚才我们得到的[latex]\dot{R}(t)=\omega*R(t)[/latex]：

也就是

注意到刚才那个*号的定义：[latex]\omega(t)*a=\omega(t)\times a[/latex].所以

这样就把刚体上质点的速度分解成了两个分量：一个是平动速度v(t)，一个是转动速度[latex]\omega\times(r_i(t)-x(t))[/latex]。注意记得x(t)是刚体质心位置。

2.6 质心

在world space下的质心就是

注意我们定义body space下质心在原点：

也就是说[latex]\Sigma m_ir_{0i}=\mathbf{0}[/latex].

当然可以算出来t时刻的质心就是x：

这里还要注意

2.7 力和力矩

记i号质点受力[latex]F_i(t)[/latex]，其力矩为

可以直觉地认为：力矩[latex]\tau_i(t)[/latex]的方向，就是刚体由于受力[latex]F_i(t)[/latex]而沿其旋转的坐标轴。

总力矩当然就是所有的加起来：

总受力F(t)不包含所有受力怎么分布的信息。但是总力矩可以。

2.8 线性动量

质点的线性动量p定义为

总的线性动量就是所有质点求和：

把每个质点的速度公式写出来就是：

根据(2-20)，最后那一项实际是零。所以其实刚体总线性动量就是刚体总质量乘以质心速度：

这样就可以用线性动量计算速度的变化率：

之所以这么写，是因为线性动量的变化率很好求，就是总受力：

这其实就是牛顿第二定律：

2.9 角动量

记住线性动量的定义是P(t)=Mv(t).我们定义总角动量L(t)为[latex]L(t)=I(t)\omega(t)[/latex]，其中I(t)是一个3*3矩阵，称之为“惯性张量”(inertia tensor).它描述了刚体的质量是怎么绕着质心分布的。这个张量取决于刚体的方向，但不取决于其位置。角动量的变化率就是总力矩：

2.10 惯性张量Inertia Tensor

定义[latex]r_i'[/latex]是i号质点离质心的距离：[latex]r_i’=r_i(t)-x(t)[/latex]。然后惯性张量取如下形式：

当然实际写起来不是求和而是积分。这个式子很麻烦，但我们实际上并不需要每次都重新计算一次I。原因是：首先注意到

我们可以把I重写为

同时注意到

把3*3单位矩阵记作黑体1（因为I已经被用了……），那么I就是：

我们注意到，旋转矩阵[latex]R(t)R(t)^T=\mathbf{1}[/latex]。

这里[latex]r_{0i}^Tr_{0i}[/latex]是标量，所以

定义

那么我们就有

这个[latex]I_{\text{body}}[/latex]是定义在body space的。所以只需要计算一次就行了。

而，I的逆就是：

2.11 刚体运动方程

刚体的状态向量为

除此之外的几个辅助物理量可以由如下方式计算出来：

其变化率为：

3 计算Y的变化率

略（这部分都是代码）

4 四元数vs旋转矩阵

用单位四元数（四个元素，被normalize到单位长度）表示旋转比3*3旋转矩阵更好。

在刚体模拟中，不用旋转矩阵的重要原因之一，是为了避免数值漂移（numerical drift）。原因是，如果我们去更新旋转矩阵：

那么随着数值误差的积累，最后会导致R不再是一个旋转矩阵。这会让最终的结果变得奇怪（skewing effect）.

用四元数能减轻这个问题。原因是旋转有3个自由度，但旋转矩阵有多达9个自由度，四元数则只有四个。四元数的数值漂移体现在，它的长度会变得不再是1，但这个normalize一下就好了。所以我们直接用单位四元数q(t)表示刚体朝向。

下面我们把四元数[latex]s+v_x\mathbf{i}+v_y\mathbf{j}+v_z\mathbf{k}[/latex]写作二元组[s,v]。乘法操作为

绕单位方向u旋转[latex]\theta[/latex]就是：

如果q1和q2代表两个旋转，那么q2q1就是二者的复合（先转q1再转q2）.刚体四元数的变化率为（证明见附录B）：

这个乘法的意思是用四元数[latex][0,\omega(t)][/latex]乘以q。注意这个式子和之前的：

很像。

（后略去代码部分）

5 例子

略

6 非穿透约束

有两种类型的接触。第一种是在接触点，双方相向运动，这是碰撞接触(colliding contact)。这就需要立刻中断求解器，修改物体的速度，然后继续求解。这是因为在求解器里，物体运动是连续的。另外一种是在接触点，一个物体静息在另一个上，这是静息接触(resting contact)，比如一个球静止在地板上。这种情况下，需要计算一个让二者互相不穿透的力，也就是接触力(contact force)，这种情况不需要中断求解器，只需要作为一个力加到求解器上就行了。

可以二分寻找确切的碰撞时间。也可以去计算碰撞时间。

7 碰撞检测

可以用到某种扫描线法，具体不详

8 碰撞接触

只考虑点/面碰撞和边/边碰撞。边/边碰撞假设两边不平行。不考虑点/点碰撞和点/边碰撞。

计算出碰撞点的相对速率（标量）：

pa,pb就是两个物体上的碰撞点。这个速率为正就分离，不用管它。为零那就是静息接触。为负就是碰撞接触。这里先讲碰撞接触的情况。我们需要给它一个impulse，让物体速度瞬间变化。

设这个impulse为J。速度改变就是：

力矩：

如果无摩擦，impulse的方向是沿着法向的：

使用减号上标代表碰撞前的状态，如：

加号上标代表碰撞后状态：

对于无摩擦情况，有经验法则：

[latex]0\leq\epsilon\leq 1[/latex].等于1就是完全弹性碰撞。等于零就是完全非弹性碰撞。

具体计算如下。首先，碰撞后速度：

把两项写开就是：

其中[latex]M_a[/latex]就是a物体的质量，[latex]I_a(t_0)[/latex]为惯性张量。代入就是：

这是一个关于j的线性函数。注意：这里作者打错了，[latex]\frac{\vec{n}(t_0)}{M_a}[/latex]一项不和[latex]r_a[/latex]做叉积。

类似对物体b：

和a同样打错了。

二者之差就是：

这个是对的。然后为了计算相对速率，把这个和法向做点积。当然法向和自己点积之后就消掉了：

根据上面那个经验公式

反解出来j：

有一个trick：如果想要“固定”一个物体，那么让它的[latex]1/mass[/latex]为零，惯性张量为3*3零矩阵即可。

9 静息接触

接触力有三个要求：第一是能避免穿透。第二是不能把两个物体黏在一起。第三是希望在物体开始分离的时候，接触力为零。比如说一个砖块静止在地上，然后你把它提起来，提起来的一瞬间，接触力应当为零。

对于接触点i，我们规定一个“接触距离”[latex]d_i(t)[/latex]，为正代表分离，为负代表穿透，为零代表接触。静息接触的情况下，[latex]d_i(t_0)=0[/latex]，接下来应该保证它非负。

对于点/面接触，可以写出一个简单的方程：

这个n当然是面的法向。对于边/边接触也类似。我们需要让这个[latex]d_i[/latex]的导数非负。这个求导就是：

因为接触所以pa=pb，第一项为零。第二项其实就是刚才碰撞接触的那个相对速率了。对于静息接触，这个d的导数应当为零（因为二者维持这种状态）。这样，d和d的导数都是零。然后再求一下二阶导：

同样，第一项为零：

这个代表，两个物体之间的相对加速度。如果大于零就是有远离的趋势，等于零维持不变。当然不能让它小于零。现在假定物体B固定，法向固定（物体B是那个面），那么法向的导数为零，第二项也可以略去了。

总之，我们得到了一个条件：

这是对接触力的一个约束，因为d的二阶导取决于接触力。

然后再看之前所说的第二，第三个要求。第二个要求是不能粘着，这就是说这个接触力得是正的（在法向上）：

第三个要求的意思是，如果d大于零，那么f就是零，这个可以写成：

然后我们需要把每对物体之间的接触距离的二阶导表述成若干接触力的线性形式。这个是因为一对静息接触的物体可以有若干个接触点。所以，这就是一个熟悉的矩阵形式了。

notes里面描述了怎么构建这个矩阵，怎么解。本文略。

2 运动方程

文中推导的是不可压缩弹性固体的方程。目标是让弹性方程看起来尽量像流体方程一样。

(q,r,s)是固体上的曲线坐标。X(q,r,s,t)是该固体质点在t时刻的位置。M(q,r,s)是质量密度。由能量密度函数E[X]定义固体材料性质。弹性力就是它的Frechet导数：

这里[latex]\wp[/latex]是小扰动的意思。[latex]\wp\mathbf{X}[/latex]的意思就是在X基础上的小扰动。类似地[latex]\wp E[/latex]就是对总能量的小扰动。

定义体积改变比例J：

弹性固体材料不可压缩也就是：

根据最小作用量原理，物体运动应当最小化：

这个L是Lagrangian.而运动过程满足固体不可压方程(2.10)和初始/终止位置：

运动过程的Lagrangian一般是动能减去势能：

因此位移的能量最小化就是：

文中说：第一项是对t做分步积分而来。第二项需要用到F的定义，也就是弹性势能E的Frechet导数。不过看上去，这个的意思好像是，[latex]\wp\mathbf{X}[/latex]是对X做小扰动，那么[latex]\wp S[/latex]就是[latex]\wp L[/latex]的积分，计算方法就是L对X求导，然后乘以小扰动[latex]\wp\mathbf{X}[/latex]（线性近似）。

定义Eularian速度场u，那么不可压缩方程就是[latex]\nabla\cdot\mathbf{u}=0[/latex].再由定义扰动v：

也就是说，u(x,t)是t时刻x处物质点的速度。v是该物质点的扰动（和静息状态相比）。

对Eularian速度u求物质导数：

进而有加速度：

v虽然看上去是速度，但其实是位移。由于固体不可压，J应该保持不变。

令

那么[latex]\partial X/\partial a[/latex]是一个3*3矩阵。而

我们现在该式两端施加扰动算子。

有恒等式：

推导如下：

这个相当于是把[latex]A_{ij}[/latex]的cofactor（也就是伴随矩阵里的i,j项）写成[latex]\partial det(A)/\partial A_{ij}[/latex].

因此用全导数写出来就是：

根据链式法则：

变换一下：

两边取trace：

因此不可压缩方程等于：

这和Eularian速度的不可压缩方程类似：

使用Dirac函数可以得到：

这个意思可以大致理解为[latex]\wp X(q,r,s,t)=v(x,t)[/latex].

然后代入刚才的最小能量条件(2.15)：

这里面质量M和弹性力密度F是Lagrangian的，但其他是Eularian的。为了把Lagrangian消掉，我们先定义Eularian量：

然后就可以把(2.32)写作：

位移v满足v(x,0)=v(x,T)=0（注：这里没太明白）和无散条件[latex]\nabla\cdot\mathbf{v}=0[/latex].

这里使用Hodge分解：

其中

下面证明w=0.设

要求[latex]\phi(0)=\phi(T)=0[/latex]，且令(0,T)之间[latex]\phi(t)>0[/latex]，那么(2.37)变为：

由于[latex]\nabla\cdot w=0[/latex]，所以p那一项被消去了，剩下：

由于[latex]\phi[/latex]总为正，所以[latex]w=0[/latex].

然后再加上一个粘性项（像牛顿流体一样），整个不可压缩弹性固体的运动方程就是：

一定要注意，这实际上是固体的方程。等于是把用狄拉克函数把Lagrangian物理量转移到Eularian网格上，从而把固体方程写成了流体的形式。

一般来说流体还会有一个质量守恒方程：

也就是[latex]D\rho/Dt=0[/latex].这个可以从(2.44),(2.45),(2.47)推出：

(2.45)两侧对t求偏导得：

而(2.45)与u点积而得：

然后首先证明一下：

证明方法是，对于任意函数[latex]\phi[/latex]：

这对任意[latex]\phi[/latex]成立，所以(2.52)成立。这样就可以从(2.50),(2.51)推出(2.49).

3 流体和固体的交互

流体的方程其实跟前面列出的固体方程差不多。但是流体不涉及弹性势能函数E，它的Lagrangian弹性力[latex]\mathbf{F}=-\wp E/\wp\mathbf{X}[/latex]等于零，对应的Eularian值f也是零。

假设流体有固定密度[latex]\rho_0[/latex]，用固体的“额外密度”[latex]\tilde{M}[/latex]代替Lagrangian密度M：

等于是密度减去它排开水的密度。这个额外密度有可能是负的，也就是固体比流体更轻。

对于浸没在流体中的表面，描述一个点只需要两个Lagrangian变量(r,s)，那么交互方程就是：

这里面，由于排开水的体积为零，所以额外密度等于密度。

这里的方程等价于那种带pressure jump的方程。证明可见Peskin and Printz(1993)，Lai and Li(2001).

4 空间离散化

（待续）

3 简介

MPM是结合PIC和FLIP而发展的一种算法。MPM对Lagrangian mesh连通性没有要求。

和PIC/FLIP类似，MPM算法在背景Eularian网格的辅助下，隐式处理自碰撞和破碎。与传统的Lagrangian方法（例如FEM算固体）Eularian方法（流体）相比，MPM好处如下：

• 和FEM一样，MPM可以由弱动量守恒形式导出，因此求解精确。
• 可以很容易在网格和粒子上施加边界条件、固壁碰撞和外力。
• 自动解决自碰撞/接触。因为粒子的运动，是由网格上未被扭曲的节点运动插值而得。
• 自动对物质进行分割/合并，因为物质是由粒子表示的。这对流体和颗粒状材料很有用。
• 只要给粒子赋不同的性质或成分模型，就能自动解决多材质和多相耦合问题。
• MPM也能用于模拟基于网格的Lagrangian力，而不失其优点（见12.2节）。因此可以同时耦合用点表示的物体，和用mesh表示的物体，只需要在网格上求解一次，并能自动处理碰撞。

严格来说MPM是Lagrangian方法，但有一个Eularian网格用来算导数。坏处是拓扑结构变化、碰撞、超弹性材料效果不佳。好处是不用单独处理碰撞和拓扑结构变化。

4 MPM在工业界的应用

《冰雪奇缘》里的雪是用MPM做的。《超能陆战队》和《疯狂动物城》用到了MPM。

5 运动学理论

每个MPM粒子实际代表一个体积微元，而不是什么质点，原子之类的东西。也可以把它看作FEM里面四面体的顶点（见第9章）。

5.1 连续物质的运动

记[latex]\mathbf{X}[/latex]为“初始位置”，[latex]\mathbf{x}[/latex]为“当前位置”。初始位置集合记为[latex]\Omega^0[/latex]，当前位置集合记为[latex]\Omega^t[/latex]，它们都是[latex]R^d[/latex]的子集。[latex]\phi[/latex]是二者之间的映射关系：

物质点 [latex]\mathbf{X}[/latex] 的速度为

加速度为

这里的速度和加速度都是以Lagrangian方式定义的，也就是说在物质点上定义的。

5.2 形变

上面说的[latex]\phi[/latex]又叫做形变映射（deformation map）。它的Jacobian可以描述一些性质，比如说弹性。这个Jacobian一般记作F：

F一般简称为形变梯度（deformation gradient），也就是一个2*2或者3*3矩阵。只有一个特殊情况：对于三维布料/薄壳，F是3*2矩阵，因为材料是二维的。可以把F看作一个映射，对于某个初始点[latex]\mathbf{X}[/latex]，它的[latex]F(\mathbf{X},t)[/latex]是对应这块物质微元在时间t的形变梯度。

写开就是：

F的行列式的含义是微元的体积变化，一般把这个行列式记作J。如果J>1意味着体积增大，J<1意味着体积减小。

J=0意味着体积变为零（退化）。J<0意味着“反转”，比如一个三角形，一个顶点越过对边。第6章会讨论“invertible elasticity”模型。

5.3 前推和拉回

物质的初始状态[latex]\Omega^0[/latex]和当前状态 [latex]\Omega^t[/latex]在[latex]\phi[/latex]下同胚/微分同胚，因为 [latex]\phi[/latex] 是光滑的。 [latex]\phi[/latex] 存在一个逆映射。所以如果一个函数定义在[latex]\Omega^0[/latex]上，就可以把它给“转移”到 [latex]\Omega^t[/latex]（这叫‘前推/push forward’），当然反之也可（‘拉回/pull back’）.例如：有一个函数[latex]G: \Omega^0\rightarrow R[/latex]，可以把它前推成[latex]g(\mathbf{x},t)=G(\phi^{-1}(\mathbf{x},t))[/latex]，而g可以拉回为[latex]G(\mathbf{X})=g(\phi(\mathbf{X},t),t)[/latex].

有时把前推函数称为Eularian（定义在[latex]\mathbf{x}[/latex]上），拉回函数称为Lagrangian（定义在[latex]\mathbf{x}[/latex]上）。

速度和加速度都是Lagrangian的：

把它前推成Eularian的：

当然可以反过来写出拉回的结果：

使用a和v的定义可以写出：

单个分量就是：

最后一项是用了爱因斯坦求和约定，代表对j=1,2,3求和。

把(8)和(10)联立，用[latex]\phi[/latex]表示Eularian速度：

联立(11)和(15)：

可以发现

注：这里求Lagrangian速度用的是偏导数，但是求Lagrangian加速度用了全导数，所以加速度多出来一项。也就是说，Lagrangian速度和Eularian速度相同，但加速度不同。Lagrangian加速度是速度的物质导数。

5.4 物质导数

物质导数：

Lagrangian加速度是速度的物质导数：

对于任意Eularian函数，物质导数的定义是：

注意它是偏导数[latex]\frac{\partial}{\partial t}F[/latex]的前推。这个F是Lagrangian的，是f的拉回。这一点很重要。

形变梯度一般认为是Lagrangian的。不过也可以把它前推为Eularian的：

这里重复下标k代表求和。推导如下：

最后一步是由对(12)求导而得。

(23)式对后面求每个粒子F的离散化更新很重要（9.4节）。

5.5 体积和面积变化

初始体积微元dV在一段时间后变成了dv=JdV，其中J=det(F)。可以用三维长度推一推得到这件事情。

所以在初始物质上积分和在一段时间后物质上积分满足关系：

为了便于理解写成：[latex]dx=JdX[/latex].

二维有向面积：

大写S是初始状态，小写s是中间状态。

考虑一个线段微元[latex]dL[/latex]（中间状态[latex]dl[/latex]），和dS（或ds）确定一个体积：

根据上面说的[latex]dv=JdV[/latex]，所以

注意F是形变梯度。在微元内dl=FdL是对的。这对任何dL都成立。：

为什么是逆转置呢？把向量点积写成向量和向量转置相乘的形式就知道了。

如此一来，在二维表面上的积分形式如下：

其中H是h的拉回。

6 超弹性

Stress是应力，strain是应变。

超弹性材料的第一类Piola-Kirchoff应力张量P可以由应变能量密度函数（strain energy density function）[latex]\Psi(F)[/latex]表示：

分量就是：

别忘了[latex]\Psi(F)[/latex]是弹性势能密度函数（是一个标量函数）。P和它大小相同。

在工程领域经常用Cauchy应力张量[latex]\sigma[/latex]，它和P的关系是：

固体材料的表现就由形变映射[latex]\phi[/latex]、Cauchy应力张量 [latex]\sigma[/latex] （当然也可以用P）这两个量确定。

对于刚体而言，形变映射就是旋转加平移：

超弹性材料可以变形，但是会抵抗变形，这由[latex]\Psi[/latex]描述。注意刚体的F=R，这是一个正交矩阵。所以对于超弹性材料，当F非正交的时候，会产生应力。正交矩阵，比如R的特性是[latex]R^TR=I[/latex]，所以我们可以写成：

这样看起来较为清楚，因为正交是[latex]F^TF=I[/latex]的特殊情况，此时[latex]\tilde{\Psi}[/latex]达到最小值.这个[latex]F^TF[/latex]经常记为C（右Cauchy-Green张量）。如果材料是各向同性的，那么可以把能量写成[latex]F^TF[/latex]特征多项式的三个系数（或者叫做各向同性不变量/isotropic invariants）：

注：各向同性的意思，是不受坐标系旋转影响。也就是说乘以一个旋转矩阵结果一样。

这个结果可以简写为：

这里面[latex]F=U\Sigma V^T[/latex]是F的“Polar SVD”.这里面U和V是旋转矩阵。而F的极分解F=RS可以通过[latex]R=UV^T[/latex]和[latex]S=V\Sigma V^T[/latex]得到 。这个R是距离F最近的旋转矩阵。S是对称阵。

各向同性不变量通常可以用奇异值表示。

6.2 新胡克模型

弹性势能密度函数：

d是维度，[latex]\mu[/latex]和[latex]\lambda[/latex]由杨氏模量E、泊松比[latex]\nu[/latex]确定：

记得J是形变梯度F的行列式。如果F是旋转矩阵，[latex]\Psi=0[/latex].如果没有反转（J>0），那么[latex]\Psi(F)\geq 0[/latex].

为了计算力，需要给出P的形式：

为了隐式积分，需要给出[latex]\frac{\partial P}{\partial F}[/latex]，这个在6.4节讲。

6.3 改进Corotated本构模型（Constitutive Model）

这个模型是基于SVD定义的。“改进”的意思是，它是在计算机图形学常用的corotated线性弹性模型制上进行改进。做polar SVD分解[latex]F=U\Sigma V^T[/latex]，那么改进corotated模型为：

这里的[latex]\sigma_i[/latex]是SVD分解里面[latex]\Sigma[/latex]矩阵的对角元。J等于各个 [latex]\sigma_i[/latex] 之积（行列式）。把右端[latex]\mu[/latex]项展开：

可以证明：

其中F=RS是极分解。给定SVD分解[latex]F=U\Sigma V^T[/latex]以后，[latex]R=UV^T[/latex]，[latex]S=V\Sigma V^T[/latex].

因此可以计算出P：

二阶导数稍微麻烦一点。

后面的太麻烦，先不管。6.4节其余部分（包括6.4.1和6.4.2）暂略。

6.5 雪的塑性（暂略）

7 物理方程

Lagrangian方程：

还记得J是形变梯度矩阵F的行列式，代表体积变化。R是Lagrangian意义下的密度。它和Eularian意义下的密度[latex]\rho[/latex]关系是：[latex]R(\mathbf{X},t)=\rho(\phi(\mathbf{X},t),t)[/latex].这个[latex]\phi[/latex]就是形变映射。注意质量守恒可以写作[latex]\frac{\partial}{\partial t}(R(\mathbf{X},t)J(\mathbf{X},t))=0[/latex].

Eularian方程：

这里v是Eularian速度，D是物质导数：[latex]\frac{D}{Dt}=\frac{\partial}{\partial t}+v\cdot\nabla^{\mathbf{x}}[/latex].

在Lagrangian方程中使用第一Piola-Kirchoff张量P更方便，在Eularian方程中则是Cauchy张量[latex]\sigma[/latex]更方便。

7.1 质量守恒

[latex]\rho(\mathbf{x},t)[/latex]是Eularian密度，R是Lagrangian密度（也就是[latex]\rho[/latex]的拉回）。二者满足关系：

Eularian密度的定义：

这个[latex]B_{\epsilon}^t[/latex]是t时刻，以[latex]\mathbf{x}[/latex]为秋心，半径为[latex]\epsilon[/latex]的球。注意这里面的体积微元实际上都是在初始状态意义上讲的。取极限就是：

注意[latex]J(\mathbf{X},0)=1[/latex]（F是单位阵）。所以这个式子可以写成：

为了推导Eularian方程，首先把它展开：

继续展开：

其中使用了行列式对矩阵求导的公式：

其中这个[latex]F^{-T}[/latex]的含义是[latex](F^{-1})^T[/latex].

(99)第一个等号，脚标的意思是求和。第二个等号，是考虑到F是形变梯度[latex]\phi[/latex]的Jacobian矩阵，也就是说[latex]F_{ij}=\partial\phi_i/\partial X_j[/latex].而Lagrangian速度[latex]V=\frac{\partial\phi}{\partial t}[/latex].所以这个[latex]\partial F_{ij}/\partial t[/latex]就是[latex]\partial\phi_i/(\partial t\partial X_j)=\partial V_i/\partial X_j[/latex].

第三个等号是基于5.2节的(10)式的，使用链式求导法则：[latex]\frac{\partial V_i}{\partial X_j}=\Sigma_k\frac{\partial v_i}{\partial x_k}\frac{\partial x_k}{\partial X_j}= \Sigma_k\frac{\partial v_i}{\partial x_k}F_{kj}[/latex].当然(99)中用爱因斯坦求和约定省去了求和号。

第四个等号，我们考虑其中固定i,k，对j求和的三项，会发现实际上是[latex]F^{-1}[/latex]的一列乘以F的一行，所以就是[latex]\delta_{ik}[/latex]，也就是i=k时为1，否则为0.

联立(97),(98),(99)：

这个J实际可以约掉。然后对两项都做前推。第一项的前推就是[latex]\rho[/latex]的物质导数，而R的前推就是[latex]\rho[/latex]，结果就是：

注意后面那项看起来很像物质导数的对流项，但其实不是。写开之后，实际上有三项。

7.2 动量守恒

作用在物质上的力可以分为体力（如重力）和面力（基于应力，stress-based）。

面力意味着有一个应力场（柯西应力张量）：

t是面力。

动量守恒公式：

这个很好理解：体积微元内总动量的变化率等于总受力。第二个等号，用到了上面的质量守恒：[latex]R(X,t)J(X,t)=R(X,0)[/latex].

(105)的最后一个等号可以写成：

(106)左端前推：

最后一个等号用了高斯公式。然后对微元取极限就是：

当然也可以对(106)右端做拉回。首先注意到：

其中[latex]F^{-T}[/latex]这些是5.5节中推出来的边界积分的结果。

因为第一类Piola Kirchoff应力张量P和Cauchy应力张量的关系是[latex]P=J\sigma F^{-T}[/latex]，所以有：

结合(106)，就得到了Lagrangian形式的动量守恒方程：

取极限：

这里面的[latex]F^{ext}[/latex]是[latex]f^{ext}[/latex]的拉回。

7.3 弱形式

首先忽略外力，动量守恒方程为：

也就是：

这里[latex]R_0=R(X,0)[/latex].然后[latex]P_{ij,j}[/latex]的意思是[latex]P_{ij}[/latex]对第j坐标轴求导，当然用爱因斯坦约定隐含了对j求和。

对于任意函数[latex]Q(\cdot,t):\Omega^0\rightarrow R^d[/latex]，把它和(113)式做点积并积分起来，得到弱形式：

注意，这里似乎应该是只对j求和，不对i求和。

其中[latex]P_{ij}N_j[/latex]应该被作为边界条件给出。设它是T(X,t)，那么对于任意函数Q，应有：

在MPM中，应力导数是在当前状态下计算的，所以可以把涉及应力的积分前推至Eularian的。令q是Q的前推，即[latex]Q(X,t)=q(\phi(X,t),t)[/latex]，反过来[latex]q(x,t)=Q(\phi^{-1}(x,t),t)[/latex]，则：

回忆起Cauchy应力张量的定义：

令t表示Eularian上边界的单位面积受力（即T的前推），则(118)变为：

我们一项一项分析。

等号左端，[latex]\rho(x,t)J(X,t)=R(X,0)[/latex]，而[latex]dx=J(X,t)dX[/latex]，合起来就是[latex]\rho(x,t)dx=R(X,0)dX[/latex].

右端第一项，本质上就是边界受力求和，是面积上加和，至于面积变化均摊的部分实际上是算到t里面了。

右端第二项就是上面说的，Cauchy应力张量。这里应该是对k求和，而不对i求和。

(120)对任意的q都成立。我们把(120)称作Eularian下的弱受力平衡形式。MPM基于这个式子做离散化。

8 物质点

在物质点上储存质量、速度、位置。但是应力是在Eularian网格上算的，所以需要一个转移的方法。

8.1 Eularian插值函数

SPH是在粒子周围定义核函数。MPM不一样，是在Eularian格点上定义核函数（插值函数）。在格点i上定义插值函数[latex]N_i(x)[/latex]。这里的i其实应该是一个向量（原文中写的是黑体，我懒得每次mathbf了）。对于某个粒子点[latex]x_p[/latex]，常简记[latex]N_i(x_p)=w_{ip}[/latex].

其中p是这个粒子，h是网格边长。之后简记：[latex]w_{ip}=N_i(x_p)[/latex]，[latex]\nabla w_{ip}=\nabla N_i(x_p)[/latex].

当然这里核函数N的选择很关键。我们喜欢用张量积样条（tensor product spline）.分段线性核函数对于FLIP很简单，但在这里不适用。原因有二（[Steffen et al., 2008]）：第一，如果用分段线性，[latex]\nabla w_{ip}[/latex]就不连续了，产生不连续的受力。第二，用分段线性，则在离格子远处，[latex]w_{ip}[/latex]本身很小，但[latex]\nabla w_{ip}[/latex]却很大，这不科学。MPM中使用的核函数最好具有C1连续性。二次B样条计算更快，三次B样条计算慢一些，但是效果更好。三次核函数为：

h是网格边长。二次核函数为：

如图3所示。

虽然理论上不稳定，但某些特定情况下也可以用线性核：

计算梯度的方法就是每一维单独算：

8.2 Eularian/Lagrangian质量

每个物质点代表一小块物质，质量就是这一块物质的质量。

把质量和动量（由此可以计算出能量）从物质点转移到Eularian网格上的方法就是：

[latex]m_i[/latex]就是坐标为i的网格的质量。为了质量守恒，核函数应当满足单位分解（unity partition），质量守恒率写出来如下：

8.3 Eularian/Lagrangian动量

类似地，可以把粒子的动量[latex]m_pv_p[/latex]转移到网格上（注：我们最后实际使用10.1节所述的APIC方法）：

类似地，它同样满足守恒律：

用动量除以质量就是速度：

这里没有爱因斯坦求和约定。

把Eularian值转移到Lagrangian值的过程类似。不过不需要转移质量，因为Lagrangian质量总是不变的。速度可以通过插值得到：

动量守恒律：

满足动量守恒律的关键一步，是由于格点质量乃从Lagrangian质量插值而得。

9 离散化

9.1 时间离散化

根据之前的弱形式：

这个[latex]t^n[/latex]是当前时间步。把加速度A离散化，左端前推：

这个能成立，要记得之前说过的，[latex]\rho dx=RdX[/latex].另外式中[latex]v^n[/latex]代表第n个时间步的速度。

为了不把代表第几维的i和代表网格坐标的i混淆，接下来我们用希腊字母[latex]\alpha,\beta,\gamma[/latex]代表第几维。由此重写(134)为：

9.2 空间离散化

类似FEM，把之前的q和v都用格点插值代入：

适当省略后简写为：

则受力平衡式(135)化为：

即：

其中质量矩阵

中间积分项的含义是：用格点i在x处的核函数，乘以格点j在x处的核函数。

拉回得到：

这是一个对称半正定矩阵，它可以写成[latex]BB^T[/latex]，其中[latex]B_{ip}=\sqrt{m_p}N_{ip}[/latex].但这个矩阵是奇异矩阵，我们使用mass lumping技术来构造一个正定对角矩阵来近似它。

(141)对所有q成立。我们取q为：

注意它是定义在格点上的。黑体i就是格点坐标。

则有：

这就是离散的受力平衡方程。这个可以作为显式Eularian更新动量的方程。

经常使用mass lumping简化。方法是把[latex]m_{ij}[/latex]的每一行加起来放在对角线上。这样就是一个对角矩阵。格子i对应的行之和[latex]\hat{m_i}[/latex]就是：

同样做简化，类似(143)：

原因是[latex]m_p\approx V_p^0R(X_p,0)[/latex].也就是说[latex]\hat{m_i}[/latex]可被自然地近似为[latex]m_i[/latex].

把[latex]m_iv_i^n[/latex]写成动量[latex](mv)_i^n[/latex]，则离散方程近似为：

左边是动量的改变，右边近似是受力。假定我们已经对每个Lagrangian点[latex]x_p^n[/latex]计算出了张量[latex]\sigma_p^n\approx\sigma(x_p^n,t^n)[/latex]，则：

其中[latex]V_p^n[/latex]是[latex]B_{\delta x,p}^{t^n}[/latex]，也就是粒子p在时间步[latex]t^n[/latex]占的体积。

用两种方法计算[latex]V_p^n[/latex]。第一，注意到

可近似估计

故

第二种方法如下。物质点p的形变梯度矩阵近似为：[latex]F_p^n\approx F(X_p,t^n)[/latex].那么其体积变化[latex]J_p^n=det(F_p^n)[/latex].假如我们知道初始体积[latex]V_p^0[/latex]，则：

使用第一类Piola Kirchoff stress：[latex]\sigma=\frac{1}{J}PF^T[/latex]，则可以把(147)和(148)中的[latex]\sigma[/latex]用P表示：

由于大部分本构模型都是指定P，所以这个写法很有用：

9.4 形变梯度更新

形变梯度较难计算。使用如下方程（见5.4节）：

注意这个式子很神奇，把形变梯度的时间导数转换成了速度的空间导数。假设只有第n+1步的Eularian速度。 则有：

做时间离散化就是：

那么更新方程就是：

当然，把速度的插值公式写出来：

结果就是：

这就是给定[latex]v_i^{n+1}[/latex]和[latex]F_p^n[/latex]后更新[latex]F_p^{n+1}[/latex]的方法。

9.5 作为势能梯度的力

对于超弹性材料，第一类Piola Kirchoff应力张量由弹性势能密度定义：

如第6章所述，某节点处的弹性势能为：

其中[latex]\hat{x}[/latex]是所有Eularian格点[latex]\hat{x}_i[/latex]拼接而成的向量。这里相当于是，临时认为Eularian网格也可以动，所以用的是Eularian格点。特别地：

那么粒子处的形变梯度更新方式为：

在[latex]\hat{x}_{i\alpha}[/latex]处对势能求导：

此外，根据上面F的方程，求导有：

代入(167)得：

和9.3节的结果对比，发现网格受力为：

所以动量更新方程就是：

10 显式时间积分

10.1 APIC转移

可以用Affine Particle-In-Cell把物理量从粒子转移至网格。忽略代表时间的上标。

转移方程：

其中[latex]B_p[/latex]是存在粒子处的一个矩阵。[latex]D_p[/latex]为：

从网格转移至粒子的方程：

使用[latex]D_p[/latex]可以轻松表示MPM常用的二次（[latex]D_p=\frac{1}{4}\Delta x^2I[/latex]）或三次（[latex]D_p=\frac{1}{3}\Delta x^2I[/latex]）插值权值。

如果使用线性样条（虽然不推荐），为了避免[latex]D_p[/latex]成为奇异矩阵，需要改为：

而从网格到粒子的转移方程就是：

此时，[latex]C_p[/latex]是[latex]x_p[/latex]处的Eularian速度梯度。我们在更新[latex]F_p[/latex]的时候已经计算过了。

10.2 形变梯度矩阵的更新

根据网格速度场[latex]v_i^{n+1}[/latex]，则形变梯度矩阵F的更新方程为：

这个就是9.4节的结果。

10.3 状态更新

显式时间更新：

根据(182)，可以把节点位置想象成关于节点速度的函数：

那么：

即：显式的力就是假设节点没有移动时的力。使用(182)，可以把(181)重写为：

这样[latex]F_p[/latex]就是节点位置的函数。

10.4 力

MPM在格点上储存力。正如之前所述，对于超弹性固体，这个力既可以用动量守恒方程的弱形式推出，也可以用弹性势能梯度推出。后者比较简单。

首先假设我们基于网格定义了超弹性势能密度。总势能为

其中[latex]V_p^0[/latex]是粒子p的体积。

节点上的弹性力等于势能的负梯度。由(186)，形变应力为：

这是基于粒子计算节点处弹性力的方法。用辛Euler方程，可以写成：

完全取决于已有的插值权值，以及粒子的物理量。如果不用能量密度，也可以用Cauchy应力张量表示这个力：

10.5 MPM整体流程

下面描述基于辛Euler方法的MPM。假设已经初始化了所有粒子（赋予质量[latex]m_p[/latex]，体积[latex]V_p^0[/latex]，初始位置[latex]x_p[/latex]，初速度[latex]v_p[/latex]，仿射矩阵[latex]B_p[/latex]和其他材料相关参数）。

1. 物理量从粒子转移至网格(Particle to Grid, P2G)。使用10.1节所述的APIC公式，把粒子上的物理量（即：质量和动量）转移至网格上。
2. 计算格点速度。[latex]v_i=\frac{m_iv_i}{m_i}[/latex]。对于质量为0的格点，把速度设为0.一般认为会在这种地方设置浮点数阈值，用来避免除以零。不过其实直接和浮点数0.0f相比也没什么问题，因为就算质量很小，导致除法运算得到一个很大的速度，但最后向粒子转移的时候，也要被这个质量再放缩回去。实际上在这里设置浮点阈值反而有害，会导致总动量不守恒之类的问题。
3. 确定自由度。这一步对于效率很关键。把质量非零的格点标记出来。其他节点没有自由度，保持不变，不在求解器中作为未知量出现。
4. 使用(189)式，显式计算格点上的力[latex]f_i^n[/latex]。
5. 使用(183) 式更新网格速度。这一不需要考虑边界条件和物体碰撞。对于显式积分而言，每个格点处的速度都可以被Dirichlet边界条件，或者刚体碰撞来单独设置。处理碰撞的细节具体见12.1节。
6. 使用(181)式更新粒子的形变梯度。注意我们并不实际把格点移到计算出来的新的位置，这只是一种假设，用来计算速度。
7. 把物理量从网格转移至粒子。这一步使用10.1节所述方法，计算出新的粒子速度[latex]v_p^{n+1}[/latex]和新的仿射矩阵[latex]B_p^{n+1}[/latex]。
8. 粒子advection。基于速度advect粒子：[latex]x_p^{n+1}=x_p^n+\Delta tv_p^{n+1}[/latex]。注意只有APIC是这么advect的。在FLIP或者混合PIC/FLIP方法中，advect方程是[latex]x_p^{n+1}=\Sigma_ix_i^{n+1}w_{ip}[/latex]。（对于PIC或者APIC，这两个式子相同）。

11 隐式积分

隐式和显式的主要区别是更新网格速度场。

11.1 力的导数

总的弹性势能可以由势能密度积分而得：

对应力做空间离散化，等于在格点上对这个势能密度做离散近似，然后求导。当然，我们并不真正改变格点的位置。设格点i的位置是[latex]x_i[/latex]，那么它下一步“应该”在的位置是[latex]\hat{x_i}=x_i+\Delta tv_i[/latex]，为简单起见，略去下标i，那么总势能就是：

这个[latex]V_p^0[/latex]就是粒子p初始时占据的体积，[latex]\hat{F}_p[/latex]为：

那么应力离散化为：

这个[latex]f_i(\hat{x})[/latex]就是格点i的弹性力。写成Cauchy应力张量[latex]\sigma_p=\frac{1}{J_p^n}\frac{\partial\Psi}{\partial F}(\hat{F_p(\hat{x})})(F_p^n)^T[/latex]：

这个[latex]V_p^n=J_p^nV_p^0[/latex]，就是时间步n时粒子p的体积。

我们使用势能的Hessian矩阵来做隐式更新。对于任意小扰动[latex]\delta u[/latex]的Hessian矩阵为：

其中

这里面冒号的含义是：A=C:D意味着[latex]A_{ij}=C_{ijkl}D_{kl}[/latex]（隐含对kl的求和）。

也可以用格点的方式推导。注意这个力在(188)中已经被推导出来，求导有：

其中[latex]\frac{\partial^2\Psi}{\partial F^2}[/latex]可以由[latex]\Psi(F)[/latex]导出。

11.2 反向Euler系统

和(183)式（[latex]v_i^{n+1}=v_i^n+\Delta tf_i(x_i^n)/m_i[/latex]）不同，反向Euler的更新方式为：

其实就是把第n步的受力换成了第n+1步的受力。也就是说，受力会隐式地依赖于网格的运动。这允许我们把时间步变得更大。

重排后得到一个非线性的格点速度公式：

11.3 牛顿法

可以用牛顿法解(200)式。初始猜测可以取上一步的速度。然后迭代步为：

每一步都是求解一个线性系统，这可以用诸如MINRES的Krylov求解器。在牛顿迭代中，[latex]F_p[/latex]需要根据速度更新。需要储存上一个时间步的[latex]F_p^n[/latex]，以备迭代使用。

在许多情况下，只需要做一步牛顿迭代，效果就非常好，能以比显式方法多出数量级的时间步做到相同效果。这几乎和所谓“半隐式”MPM积分效果相同。

11.4 线性力

弹塑性响应可以由网格位置的虚拟“更新”确定。把这种虚拟更新写成[latex]\hat{x}[/latex]关于v的函数：[latex]\hat{x}=\hat{x}(v)[/latex]。我们使用如下记号：[latex]f_i^n=f_i(\hat{x}(0))[/latex]，[latex]f_i^{n+1}=f_i(\hat{x}(v_{n+1}))[/latex]，[latex]\frac{\partial^2 e^n}{\partial\hat{x}_i\hat{x}_j}=-\frac{\partial f_i^n}{\partial \hat{x}_j}=-\frac{\partial f_i}{\partial \hat{x}_j}(\hat{x}(0))[/latex].

由此写出隐式速度更新：[latex]v_i^{n+1}=v_i^n+\Delta tm_i^{-1}((1-\beta)f_i^n+\beta f_i^{n+1})=\approx v_i^n+\Delta tm_i^{-1}(f_i^n+\beta\Delta t\Sigma_j\frac{\partial f_i^n}{\partial \hat{x}_j}v_j^{n+1})[/latex].这是一个求解[latex]v_i^{n+1}[/latex] 的对称方程组 ：

其右端就是：

参数[latex]\beta[/latex]取0时就是显式，取0.5时就是梯形公式，取1就是反向欧拉。这相当于以全0场为初值，做一次牛顿迭代。因此，这样每步只需要解一个线性系统。MPM半隐式时间步见图4。

11.5 基于优化的时间积分

Newton-Raphson求解器性能比显式时间步好得多，但仍然需要较小的时间步才能达到稳定。

如果力是由势能函数得到的（这在大多时候成立），那么就可以让牛顿法在和CFL条件相近的任意[latex]\Delta t[/latex]内收敛，方法是把远方程写成一个优化方程。

实际上，求解：

等价于最小化：

这允许我们采取大得多的时间步，节省大量运算资源。具体情况见[Gast et al., 2015].

12 其他

12.1 物体碰撞

在使用受力更新网格速度后，我们就处理碰撞问题。对于半隐式积分，这涉及到修改线性方程组的右端，并且消除涉及到碰撞的格点的自由度。可以在更新粒子坐标之前，对粒子速度[latex]v_p^{n+1}[/latex]做一次碰撞检测，用来解决插值引入的新误差。所有碰撞都是非弹性的。注意基于粒子的碰撞会让粒子上的形变梯度张量自相矛盾。因此在模拟中，只应该在必要情况下使用这种碰撞。

我们使用level set表示碰撞的物体。那么碰撞检测就很简单（[latex]\phi\leq 0[/latex]）。首先计算出相对速度[latex]v_{rel}=v-v_{co}[/latex]，这个co是碰撞物体的意思。如果两个物体互相分离，那么不加碰撞。记[latex]v_t=v_{rel}-nv_n[/latex]为相对速度的切向分量。如果需要一个“sticking impulse”（[latex]||v_t||\leq-\mu v_n[/latex]），那么直接设[latex]v’_{rel}=0[/latex]，这个撇代表已经处理过碰撞了。否则就用使用动态摩擦，更新为[latex]v’_{rel}=v_t+\mu v_nv_t/||v_t||[/latex]，其中[latex]\mu[/latex]是摩擦系数。最后基于相对速度去更新速度：[latex]v’=v’_{rel}+v_{co}[/latex].

我们有两种碰撞物体：刚体和软体。刚体就储存一个固定的level set，和刚体的位置参数，用来计算每个点的[latex]\phi[/latex]，法向和物体速度。对于软体，在关键帧加载level set，并使用[latex]\phi(x,t+\gamma\Delta t)=(1-\gamma)\phi(x-\gamma\Delta tv_{co},t)+\gamma\phi(x+(1-\gamma)\Delta tv_{co},t+\Delta t)[/latex]做时间插值。但速度则用[latex]v_{co}=(1-\gamma)v(x,t)+\gamma v(x,t+\Delta t)[/latex]，忽略位置变化。

最后考虑粘性碰撞。此时Coulomb摩擦是不够的，因为法向相对速度应该是0（垂直）或者正数（悬挂在下面，然后由于重力又分开）。这里我们无条件地设置[latex]v’_{rel}=0[/latex]。也就是说，这等价于Dirichlet边界条件。

12.2 Lagrangian力

在传统MPM模拟中，使用(188)式计算力。在粒子处储存形变梯度有一定好处，例如不用考虑拓扑。但是在用MPM模拟传统的软体，比如橡胶之类的东西时，使用基于mesh的FEM求解器更有道理，因为它提供了连通性信息，能更精确地计算形变梯度矩阵，并且有利渲染。在此情况下，MPM可以和基于mesh的方法轻松整合在一起。

我们可以使用任何Lagrangian受力模型（例如弹簧，有限元等等），只要能写出每个点处的总势能[latex]W(x_p)[/latex]。基于这个模型计算出受力[latex]f_p=-\frac{\partial W}{\partial x_p}[/latex]。设对于任何粒子处的位移[latex]\delta u_q[/latex]，都可以计算相应受力[latex]\delta f_p=\Sigma q\frac{\partial f_p}{\partial x_q}\delta u_q[/latex].这些是纯Lagrangian的，可以用常规方法做，我们就不详述了。

然后考虑如何把它们加到网格上。我们需要根据粒子处的受力[latex]f_p[/latex]计算处网格受力[latex]f_i[/latex]，并根据[latex]\delta u_i[/latex]计算[latex]\delta f_i[/latex]。考虑之前粒子-网格之间互相转移的方法，我们得到[latex]x_p=\Sigma i w_{ip}^nx_i[/latex].使用链式法则：

这里有三个隐含求和，不过计算起来并不难，可以一个一个做。这些是网格上的受力，那么MPM和Lagrangian方法就可以整合起来。把每个粒子标记为MPM粒子或者是mesh粒子。我们不会使用储存在mesh粒子上的形变梯度矩阵，因为直接在mesh上计算即可。这就能把Lagrangian方法的精确表面追踪和Eularian网格的自动处理碰撞给结合在了一起。

参考文献（略）

摘要

本文提出一种流固耦合算法，能处理用三角面片表示的，无穷薄的固体。经典的流固耦合算法在三维格子上表示固体，但薄片没有内部区域，这就需要新方法。我们使用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步固体的位置，计算出它的等效速度。使用这个等效速度，对流体的中间速度场做投影，使之满足无散条件。

然后重复上面这个过程。

7 例子（略）

8 结论和未来工作（略）

9 致谢（略）

10 参考文献（略）

单向耦合

如果固体足够小，就可以忽略固体对流体的作用力，只考虑流体对固体的作用力。

固体受到的力分为两种：压力和粘性力。在固体非常小的时候，粘性力占主导，某点处的粘性力等于粘性张量乘以法向，总的受力就是：

这里S代表固体的体积，法向指向固体外部。由于是单向耦合，因此不在固体表面应用无滑移边界条件，也就不能直接把流体的粘性张量拿来用。那怎么办呢？理论上讲，固体周围应该有边界层，应该拿边界层内的速度场算粘性张量，但边界层的流场问题还没有解决。因此书中就搞了一个简单的式子：

这里面F_i是i号固体 （非常小，近似粒子） 受到的粘性力/拖曳力，u是流体的速度场，u_i是i号固体的速度。D应当和流体的动力粘度成正比，可能根据每个固体的直径、截面积而改变。有时候甚至可以搞个非线性的，让D和u-u_i的模长成正比。

对于薄片状固体，有可能会约束速度的法向分量和固体相同，只在切向上加一个很小（甚至可能是零）的粘性力。

粘性力对固体产生的总力矩是：

其中x_i是i号固体的重心。和前面类似，没有确切公式，同样给一个简单的：

其中[latex]\Omega[/latex]是固体的角速度，[latex]\omega[/latex]是流体的旋度。E也是自己调，它可以是个矩阵。

压力的形式则简单一些：

如果固体足够小，直接在其重心处计算[latex]\nabla p[/latex]，然后乘以体积就可以了。

压力产生的总力矩是：

如果p足够平滑（接近于常数或者线性函数），这个积分的结果就可以忽略。这里注意：对于静水，压强在水面处并不平滑，因为水面上方空气都是0，但水下就是线性的，所以一半浮在水面上的物体会受到压力矩。对于这种情况，就需要对水下部分做积分（或者估算）了。

弱耦合

如果固体和流体能互相影响，那么就需要设法同时模拟二者。一种常见的做法叫做弱耦合（weak coupling）。这种方法是：交替进行流体步和固体步。最简单的算法如下：

1、流体对流步，更新固体位置和方向。
2、把非耦合的力加在流体和固体的速度场上（例如重力，内部弹性力）。
3、求解流体的压力方程，并据此修正流体速度，使之满足无散条件，固定固体速度，在其边界上对流体施加固壁边界条件。
4、根据上面求得的流体压强、碰撞、接触等，更新固体速度。

当然刻意弄得更复杂，比如多交替几次流体/固体步。但是弱耦合的本质不变：算流体步（求解压力方程）的时候假定固体速度不变，算固体步的时候假定流体压强不变。

其中，固定固体算流体就相当于流体加边界条件，在此不赘述。而固定流体算固体，实际就是计算固体受到的力和力矩，如之前单向耦合中所述。只要mesh足够细，例如和[latex]\Delta x[/latex]在同一量级上，那么可以直接做数值积分，比如说把三角形面片中间的压强乘以面积。粘性力也可以类似这么做。

但在一些情况下，数值积分可能不好操作，比如绳子这样的细长物体。这时候，可以把面积分转化成体积分，例如总压力：

总力矩：

其中[latex]x_C[/latex]是重心。体积分可以由若干网格上的值加权求和得到。

有若干文章，如：[Takahashi et al. 02, Guendelman et al. 05]采用了这种方法。这种做法在软件工程上有优势（模块分离），但是有一些问题。例如，如果初始状态是一个固体一动不动地漂浮在水面上，那么固体的速度会在第2步中增加[latex]\Delta t\vec{g}[/latex]，随后在第3步中，求解器会把这个速度作为边界条件，那么在完成速度校正后，流体的速度将非零。最终，压强场并未抵消固体的速度，固体会下沉一点，流体开始运动。这样的误差和时间步大小成正比，减小时间步可以降低误差，但计算开销会变大。

浸没边界法（Immersed Boundary Method）

浸没边界法是一个稍强一点的算法（可以看经典文章[Peskin 02]）。这个算法赋予流体压力求解器一些“余量”，来改变固体速度。方法是暂时假装固体也是流体（但密度不同）。它不是简单地把固体速度作为求解流体压强的边界条件，而是把固体的质量和速度加到流体格子上，并在整个计算域内求解压强。每个节点处都有流体和固体的权值，加权确定平均密度和平均速度（这种做法和[Guendelman et al.05]类似）。

一种与此相关的方法是[Carlson et al. 04]提出的，他们在求解压强时，假设固体的密度和流体相等，然后在另一步中单独添加浮力。该文章没有在固体表面上施加压力，而是直接把固体占据的所有格子的速度平均起来，得到固体的速度。如果流体和固体的密度差别不太大，这种方法效果就很好。

对于无粘流体，把固体和流体速度做平均，可能会产生数值耗散，因为在固体边界上为滑移边界条件，流体固体速度可以不同。而做平均相当于无滑移边界条件，让流体和固体速度相等。在这种情况下，建议把流体的切向速度外推（extrapolate）到固体格子里，只对固体的法向速度做平均。对于非常薄的固体（例如布料），这种做法就变得很简单，因为固体太薄就不需要外推了，把法向/切向速度分开即可。

相比前面的弱耦合法，浸没边界法能解决一些问题，但仍然不能处理所有情况。例如，在模拟一个静止的，漂浮在水中的固体时，结果仍然很奇怪，因为在流体的压力求解器看来，固体是流体表面一个形状奇特的波。

一般稀疏矩阵

在讲解强耦合（同时借固体和流体）之前，我们先简单介绍稀疏矩阵，在之前章节的基础上进行扩展，因为我们先前所用的稀疏矩阵储存格式不适用于固体。

储存稀疏矩阵有几种可选的数据结构。我们关注的是行压缩（compressed sparse row/CSR）格式。在CSR里，每一行都用一个数组表示，数组中包含这一行的全部非零值，和它们对应的列编号。我们使用两个版本的CSR，其一是动态版本（这样构建矩阵时，能迅速添加非零元素），其二是静态版本（做PCG效果更好）。

动态CSR里面，每一行的数组分开储存，并且记录其长度。为了支持动态插入，我们可以额外申请一些空间，在这些空间用尽后，再申请两倍大小的空间，这和C++ STL vector用的方法类似。另外，人们通常把每个数组内部按照列坐标排序，这样查询某项和把两个稀疏矩阵相加都快得多。

但是PCG的核心操作是系数矩阵乘以稠密向量，如果拿动态CSR做乘法，效率会变低，因为每一行可能会分散在不同的内存块中，导致缓存命中率极差。因此，在使用动态CSR建好矩阵后，我们应当把它转换成静态CSR。在静态CSR里只有三个矩阵，保证所有非零元素在内存中连续：

• 一个浮点数组value，包含一列一列排序后的所有非零元素。
• 一个整型数组colindex，和value长度相同，包含value中相应元素的列坐标。
• 一个整型数组rowstart，长度为n+1（假设矩阵大小为n*n），描述在value或colindex矩阵中，每一列从哪里开始。这个数组的最后一个元素是整个矩阵中的非零元素个数，它代表了最后一行的结尾。

把动态CSR转为静态CSR需要花费一定时间，不过这段时间一般是值得的，因为它大大优化了PCG中矩阵乘向量的效率，使之接近先前提到的，基于网格的储存法。

我们之前使用的，基于网格储存矩阵的方法，针对对称矩阵进行优化，只储存上/下三角阵，从而省却了一半的空间。不过这个优化对于CSR并不值得，因为它让矩阵操作变得极为复杂。

静态CSR乘稠密向量的操作很简单，见图11.1.但我们还需要把不完全Cholesky预处理子（incomplete Cholesky preconditioner）和相应的三角阵求解搬到CSR上。一般而言，先前章节中做的假设（只需要计算对角元素，然后重用非对角元素）并不成立。此外，改进版不完全Cholesky分解不能使用先前所述的循环顺序。在CSR中，计算[latex]R=L^T[/latex]更自然一些。

总之，我们需要根据如下条件，确定标准版Cholesky分解的不完全因子（regular incomplete factor）：

• R是上三角矩阵，若[latex]A_{i,j}=0[/latex]则应同时[latex]R_{i,j}=0[/latex]
• 若 [latex]A_{i,j}\neq 0[/latex] ，则应有[latex](R^T R)_{i,j}=A_{i,j}[/latex].

至于改进版Cholesky分解的不完全因子，则应满足如下条件：

• R是上三角矩阵，若[latex]A_{i,j}=0[/latex]则应同时[latex]R_{i,j}=0[/latex]
• 若[latex]A_{i,j}\neq 0[/latex]且i<j，则应有[latex](R^T R)_{i,j}=A_{i,j}[/latex].
• 每一行的和[latex]\Sigma_j (R^T R)_{i,j}[/latex]和A每一行的和[latex]\Sigma_j A_{i,j}[/latex]相等。

我们不深入解释如何使用CSR格式。不过，图11.2是构造R的伪代码（参数和之前相同），图11.3演示了如何通过求解R和[latex]R^T[/latex]，使用预处理子。

强耦合

对于刚体和无粘流体的耦合问题，强耦合（strong coupling）是目前被研究得最透彻的方法。这是本章介绍的最后一种方法。接下来我们先给出连续形式的方程，然后再把它们离散化。

首先我们对刚体定义一些符号：

• [latex]\vec{X}[/latex]是刚体重心。
• [latex]\vec{V}[/latex]是刚体平动速度。
• [latex]\vec{\Omega}[/latex]是角速度。
• [latex]\vec{L}[/latex]是角动量。
• [latex]m[/latex]是质量。
• [latex]I[/latex]是转动惯量。

（可以把它轻易扩展到多个刚体的情况）。根据先前讨论的，流体压力赋予刚体的总受力和总力矩，刚体的更新方程如下：

刚体上某点[latex]\vec{x}[/latex]的速度为：

在求解压力方程的时候，它可以作为边界条件：

另外，在流体-气体的边界，应该施加自由边界条件p=0.

此处补充一点。在边界上，把[latex](\vec{V},\vec{\Omega})[/latex]映射至法向速度场的线性算子，等于在边界上把压强映射至固体总受力和总力矩的线性算子的伴随算子（对于矩阵而言，就是转置）。从上边的式子中看去，这不很显然，但可以被证明。这个性质决定了，在离散化之后，我们的压力方程会是一个对称矩阵。

[Batty et al.07]文中给出了一个更简单的方程组。我们使用方程(11.4)和(11.5)更新固体，但并不使用上文中说的边界条件。为了简化问题，我们以最小能量原则建立压力方程：压强场应当使得动能最小化。注意刚体的动能是[latex]\frac{1}{2}m||\vec{V}||^2+\frac{1}{2}\vec{\Omega}^TI\vec{\Omega}[/latex].

接下来我们来离散化这个变分方程。我们在之前的弱耦合一节中已经讨论过，如何计算固体受压力的影响（总受力和总力矩）。下面我们只需要明确地构建一个稀疏矩阵J，使得用J乘以包含全部格点压强值的向量，就能得到受力和力矩。根据方程(11.2)，我们可以得到对应总受力的前三行。例如，总受力的x分量，[latex]F_1[/latex]如下：

其中，[latex]V_{i+1/2,j,k}[/latex]是u-格子[latex](i+1/2,j,k)[/latex]中固体的体积。因此我们有

类似地，对应总受力y和z分量的后两行分别为：

根据方程(11.3)，我们可以得到J中对应总力矩的后三行。x分量[latex]T_1[/latex]的连续形式如下：

其中，固体的重心是[latex]\vec{X}=(X,Y,Z)[/latex].

上式可以简化为

离散近似为：

其中[latex]x_i,y_j,z_k[/latex]为[latex](i,j,k)[/latex]格中心的坐标。由此可以推导出J的第四行（力矩的x分量）：

以及后两行（力矩的y，z分量）：

在固体边界之外，相邻格子的V值之差都是零，所以J是一个稀疏矩阵。它只有在固体边界附近才有非零列。

为了保持一致性，可以使用同样的V值计算固体的转动惯量，这样一来，完全浸没水中的固体就能完美地保持静止漂浮。但是，除此特定情况，这种做法大概没什么必要，我们把具体推导作为练习，留给有兴趣的读者。

为方便起见，我们把固体的平动速度[latex]\vec{V}[/latex]和角速度[latex]\vec{\Omega}[/latex]拼成一个六维向量[latex]\vec{U}[/latex].同样，把质量和转动惯量拼成一个6*6的矩阵M：

那么刚体的动能就是[latex]\frac{1}{2}\vec{U}^TM\vec{U}[/latex]，压力修正步是[latex]\vec{U}^{n+1}=\vec{U}+\Delta tM^{-1}Jp[/latex].

最后，我们来讨论离散最小化。把动能对压强求导得到：

把它和第四章末推导出的流体动能合起来，进行恰当的放缩使量纲一致（除以[latex]\Delta t[/latex]，用真实体积值代替体积分数），最终的结果是左端矩阵多了[latex]\Delta tJ^TM^{-1}J[/latex]，右端多了[latex]-J^T\vec{U}[/latex].

这正是我们需要处理一般系数矩阵的原因：[latex]\Delta tJ^TM^{-1}J[/latex]并不对应网格结构。实际上，它在边界附近格子的对应列处有着稠密子矩阵。如果你对线性代数感兴趣，可以发现它的秩至多为6，并且一定是对称半正定矩阵——所以PCG仍然能用。但是稠密子矩阵的性质可能会导致求解性能下降，可以在分解后插入一些额外项来避免这点。PCG涉及到矩阵乘稠密向量的运算，这个操作可以通过把[latex]J^TM^{-1}Js[/latex]换成[latex]J^T(M^{-1}(J^Ts))[/latex]来加速。由于秩不超过6，因此计算预处理子代价不大。

大块固体可能会和许多格子重叠，这是一个值得考虑的问题：需要用大量内存储存它们，并且PCG可能会变得很慢，由于稠密子矩阵太大。一种可能的方案是在压力求解中，把新的固体速度[latex]\vec{U}^{n+1}[/latex]作为辅助变量，得到一个稍大但更稀疏的矩阵：

这仍然是一个对称矩阵，但它可能是不定的，因此不能使用PCG。不过，对这种矩阵的求解也是一个经典问题，它有时会被叫做“鞍点”(saddle-point)矩阵（实际上，除去压力为常数的零空间，它应该是一个对称拟正定/symmetric quasi-definite矩阵），有专门的求解算法。例如，可以采用MINRES这样的迭代方法，结合形如[latex]LDL^T[/latex]（L是单位对角阵/unit diagonal，D是对角阵，其中压强值的对应项为正，刚体的对应项为负）的不完全Cholesky预处理子。如果想了解求解此类矩阵的方法，可以去看[Benzi et al. 05]的综述文章。

另一个有趣的研究方向，是考虑如何把这种方法用于关节状固体的强耦合：例如，拉格朗日乘子法可以表述为最小化动能。摩擦力也可以这么做，但它的约束是不等式，这让最小化变得很复杂。

如果需要模拟软体，上面说的最小化能量算法就不够用了。在软体-流体耦合问题中，强耦合意味着把软体本身的隐式时间步和流体的压力求解结合起来。但是，隐式求解弹性力（最小化弹性势能）不能被简单地描述为，求解一个力，来最小化能量。在写作本书时，图形学领域只有[Chentanez et al.06]研究过这个问题。他们使用非结构化的四面体网格来离散化流体，在软体边界处网格的形状和软体一致。目前，还没有使用结构化网格和固体体积分数的方法。本章至此结束。

在初始时间[latex]t_0[/latex]，在界面上放[latex]N_p[/latex]个粒子，用[latex]\mathbf{x}_p[/latex]表示粒子的位置。

显然，粒子跟着流场走，也就是：

由于是粒子，所以就是直接微分，没有什么物质导数之类的。

我们用有符号距离函数[latex]\phi[/latex] 做水平集的函数（level function）。那么根据这个[latex]\phi[/latex]就能求法向：

上面求出了粒子的时间演化方程，利用[latex]\phi[/latex]的性质可以把法向场的物质导数展开：

然后文章又对这个法向场的物质导数做了一下变换：

这个记号比较难懂，单独注明一下。首先是用了偏导的角标简写。在角标中逗号后面的做偏导。比如说[latex]\phi_{,i}[/latex]就是[latex]\phi[/latex]对第i维坐标（x,y,z之一）做偏导，如果 [latex]\phi_{,it}[/latex]就是对第i维和时间求二阶导， [latex]\phi_{j,i}[/latex]就是[latex]\frac{\partial\phi_j}{\partial\phi_i}[/latex]的意思。另外一个就是爱因斯坦求和，所有重复的角标都表示求和。比如说[latex]\phi_{,k}\phi_{,kt}[/latex]中这个重复的k就暗示一个求和，也就是说它实际是[latex]\Sigma_k\phi_{,k}\phi_{,kt}[/latex].

这个推导左端是i，因为法向场的物质导数是一个矢量，这是其中的第i维。前几个等号都可以通过微分运算法则推出来。倒数第二个等号，是专门凑了一下[latex]\phi_{,t}+\mathbf{u}_j\phi_{,j}[/latex]这个形式，把它用全微分展开即得正确性。为什么要凑这个形式，是因为文章中的水平集函数[latex]\phi[/latex]满足守恒方程：

这样凑出来的项都是零了，就只剩下两项。然后根据前面说的， [latex]\phi[/latex] 和法向的关系，很容易验证最后的结果，也就是：

[latex]-u_{j,i}n_j+u_{j,k}n_jn_kn_i[/latex]

注意这只是一维，三维加起来的结果是：

值得高兴的是，在这个式子中 [latex]\phi[/latex] 消失了。也就是说，法向的演化仅仅和速度场有关。文章采用一个三阶Runge-Kutta算法去演化法向场。

文章中采用的是有向粒子法，也就是说粒子上面携带位置信息和当地法向信息。

每一步法向都会重新normalize，放缩到长度为1.

虽然粒子提供了表面信息，但并不知道粒子之间的相邻属性，所以不能直接用于生成表面。为了生成表面，还是需要搞一个水平集函数[latex]\phi[/latex].

方法是：遍历所有粒子，确定interface cells（至少包含一个粒子的格子）。在遍历过程中，计算离每个格点最近的粒子距离。然后用这个距离和该粒子的法向做点积，即得这些点处水平集函数（有向距离）。

然后水平集函数满足Eikonan方程：

就可以把它延拓到整个计算域上。

随机采样算法

文章在此处先讨论了“表面上的两个点”之间的信息，或者说从一个点去找“下一个点”。

如图。P是一个点，然后分析表面上的另一个点的信息。这个点就是右半边图中，P右下角处的那个点。它的位置就是[latex]\mathbf{x}_P+\beta\mathbf{t}+\alpha\mathbf{n}[/latex].其中[latex]\mathbf{n}[/latex]很好理解，就是P点处的法向。至于[latex]\mathbf{t}[/latex]，则是沿表面的一个切向量。当然这样的切向量有许多，文章应该是先假设取了一个平面（图中淡黄色平面，这个平面包含[latex]\mathbf{n}[/latex]），然后全部在这个平面上讨论。

然后文章基于这个关系推导了一下。

第一步，是基于水平集函数[latex]\phi[/latex]进行推导。

把[latex]\phi[/latex]关于P点做二阶泰勒展开：

文章在这里推了一下中间那个二阶导项的具体表示：

注意这里用到了[latex]\phi[/latex]和法向的关系，尽量使用[latex]\mathbf{n}[/latex]表示。

然后把 [latex]\mathbf{x}_P+\beta\mathbf{t}+\alpha\mathbf{n}[/latex] 这个点的位置信息代入，一番推导后可得：

此处注意，新的点也应该在表面上，其[latex]\phi=0[/latex]，所以这个式子等于零。

第二步，基于几何进行推导。

还是看上面那张图的右半部分。两点之间的弧长为

这里边[latex]\tilde{\kappa}[/latex]就是和[latex]\mathbf{t}[/latex]相切这条曲线的曲率了。

进而，使用这个式子，可以把两点之间的弦长表示为

可以看到，[latex]\beta[/latex]是直角三角形的直角边，可用三角函数求其长度：

这里经过了一步三角恒等式代换和一步泰勒近似。我们可以发现这个式子给出了[latex]\beta[/latex]和s的一个关系，也就是说沿着曲面走s，能让[latex]\beta[/latex]变化多大：

然后[latex]\alpha[/latex]也可以用三角函数推导出来：

这里可以发现，[latex]\alpha[/latex]是 [latex]\alpha^2[/latex]这个量级的。把这个关系代入第一步推出来的[latex]\phi[/latex]。取一个比较小的[latex]\beta[/latex]：

然后就可以在新的[latex]\phi[/latex]表达式中消掉[latex]\beta[/latex]的三阶项：

中间这个[latex]\phi_{,ij}\mathbf{t}_i\mathbf{t}_j[/latex]也可以用法向表示。文章里说是根据第一步里面[latex]\phi_{,ij}[/latex]得来，不过我感觉好像实际的推导应该是根据[latex]\phi[/latex]和法向的关系：

首先[latex]\mathbf{n}=\nabla\phi/|\nabla\phi|[/latex]，取一维也就是说[latex]n_i |\nabla\phi|=\phi_{,i}[/latex]，两边再对第j维求导得[latex]\phi_{,ij}=n_{i,j}|\nabla\phi|+n_i\partial|\nabla\phi|/\partial x_j[/latex].而水平集函数的性质，要求[latex] |\nabla\phi| =1[/latex]，偏导总是零，所以第二项就没了，只剩 [latex]\phi_{,ij}=n_{i,j}|\nabla\phi| [/latex].

这只是一个式子，乘上[latex]t_i t_j[/latex]，再求和（别忘了上面那些都是自带爱因斯坦求和的），就得到了：

再代回去，结果就是

这样就没有[latex]\phi[/latex]了，直接用P点处的几何信息，就能用[latex]\beta[/latex]得到[latex]\alpha[/latex].

接下来，使用Frenet第一公式可以得到切向和法向的关系（公式就长这样）：

两端点积[latex]\mathbf{n}[/latex]就可以得到曲率：

第二个等号用到了Frenet第二公式。第三个等号，把ds写开推一推可以得到。

再回忆一下之前说的，新的点的位置：

[latex]\mathbf{x}_P+\beta\mathbf{t}+\alpha\mathbf{n}[/latex].

它的第i维就是：

[latex]x_i^{new}=x_i+\beta t_i+\alpha n_i[/latex].

我们刚才推出来了[latex]\alpha[/latex]的公式，略去小量，代入得

[latex]x_i^{new}=x_i+\beta t_i-\frac{1}{2}n_{k,l}t_k t_l\beta^2 n_i[/latex].

最后一项可以用曲率表示：

[latex]x_i^{new}=x_i+\beta t_i+\frac{1}{2}\tilde{\kappa}n_i \beta^2 [/latex].

注意此处文章里的公式最后一项没有[latex]\beta^2[/latex]，我认为是作者搞错了，否则说不通。

相应的，新的点处的法向是：

曲率张量是：

还记得我们之前说的，文章假定了切向[latex]\mathbf{t}[/latex]，但之前没说明怎么选取这个切向。选取的方法包括：

1、计算shape operator的特征向量，用来确定主曲率（principal curvatures）方向，并把切线设成沿这些方向。

2、简单地选定两个垂直的方向。

文章里用的是第二种。对于每一个初始种子，在其四周用上面讨论的方法，取四个采样点，也就是两个垂直的方向，及其相反方向。

根据文献[chen2008]，使用Darcy-Brinkman-Forchheimer扩展模型描述孔隙流动的方程如下：

不可压方程：[latex]\nabla\cdot\vec{u}=0[/latex]

动量方程：[latex]\rho\frac{\partial \vec{u}}{\partial t}+\nabla\cdot\left(\frac{\rho\vec{u}\vec{u}}{\epsilon}\right)=-\nabla(\epsilon p^*)+\mu\nabla^2\vec{u}-\frac{\mu\epsilon}{K}\vec{u}-\frac{\rho\epsilon C_F|\vec{u}|}{\sqrt{K}}\vec{u}[/latex]

在这里的[latex]\vec{u}[/latex]是Darcy velocity，或local average velocity，或superficial velocity。其意义是，单位截面积孔隙介质中流出的流体通量。这里的“单位截面积”是固体和液体一起算，所以这个Darcy velocity会比孔隙中流体质点的“真实速度”小一些。二者满足关系：Darcy velocity=真实速度*孔隙度（孔隙度一般记作[latex]\epsilon[/latex]）。而p*这个星号的含义是它是intrinsic也就是流体质点上的真实压力。另外，Forchheimer系数由[latex]C_F=1.75/\sqrt{150\epsilon^3}[/latex]定义，它和孔隙度一样是一个无量纲数。

我们定义如下特征量：特征速度U（定义为远处来流速度），特征尺度D（依问题不同而异，比如说管流的D可取管径，圆柱绕流的D可取圆柱直径等等）。

在这样的流动中，比较关键的无量纲参数有雷诺数[latex]Re=U\rho D/\mu[/latex]和达西数（Darcy number）[latex]Da=K/D^2[/latex]。这里我们认为流体密度[latex]\rho[/latex]是一个常量。

定义无量纲速度[latex]\vec{u}’=\vec{u}/U[/latex]，无量纲压力[latex]p’=p^*/(\rho U^2)[/latex].

不可压方程转换为无量纲形式之后就是乘了一个常数，形式不变：

[latex]\nabla\cdot\vec{u}’=0[/latex]

动量方程比较麻烦。首先无脑转换：

[latex]U\rho\frac{\partial \vec{u}’}{\partial t}+U^2\nabla\cdot\left(\frac{\rho\vec{u}’\vec{u}’}{\epsilon}\right)=-(\rho U^2)\nabla(\epsilon p’)+\mu U\nabla^2\vec{u}’-\frac{U\mu\epsilon}{K}\vec{u}’-\frac{\rho U^2\epsilon C_F|\vec{u}’|}{\sqrt{K}}\vec{u}'[/latex]

然后把每一项除以[latex]\rho U^2[/latex] ：

[latex]\frac{1}{U}\frac{\partial \vec{u}’}{\partial t}+\nabla\cdot\left(\frac{\vec{u}’\vec{u}’}{\epsilon}\right)=-\nabla(\epsilon p’)+\frac{\mu}{\rho U}\nabla^2\vec{u}’-\frac{\mu\epsilon}{\rho UK}\vec{u}’-\frac{\epsilon C_F|\vec{u}’|}{\sqrt{K}}\vec{u}'[/latex]

为了凑雷诺数和达西数，把每一项乘以D：

[latex]\frac{D}{U}\frac{\partial \vec{u}’}{\partial t}+D\nabla\cdot\left(\frac{\vec{u}’\vec{u}’}{\epsilon}\right)=-D\nabla(\epsilon p’)+\frac{D^2\mu}{\rho UD}\nabla^2\vec{u}’-\frac{\mu\epsilon D^2}{\rho UDK}\vec{u}’-\frac{D\epsilon C_F|\vec{u}’|}{\sqrt{K}}\vec{u}'[/latex]

也就是：

[latex]\frac{D}{U}\frac{\partial \vec{u}’}{\partial t}+D\nabla\cdot\left(\frac{\vec{u}’\vec{u}’}{\epsilon}\right)=-D\nabla(\epsilon p’)+\frac{D^2}{Re}\nabla^2\vec{u}’-\frac{1}{Re}\frac{\epsilon}{Da}\vec{u}’-\frac{\epsilon C_F|\vec{u}’|}{\sqrt{Da}}\vec{u}'[/latex]

至此无量纲化还未完成，因为时间、空间、U、D都是带量纲的。为了把它们干掉，按照[nithiarasu2002]文中方法，把坐标和时间也无量纲化为：[latex]x’=x/D[/latex]（y,z无量纲化方法相同，略去），[latex]t’=tU/D[/latex]。这里实际上是认为特征尺度/特征速度=特征时间（=D/U）。

如此一来，空间偏导和时间偏导也就相应地无量纲化了。特别地，无量纲的nabla算子满足[latex]\nabla’=D\nabla[/latex].如此一来我们就得到了最终的无量纲动量方程：

[latex]\frac{\partial \vec{u}’}{\partial t’}+\nabla’\cdot\left(\frac{\vec{u}’\vec{u}’}{\epsilon}\right)=-\nabla'(\epsilon p’)+\frac{1}{Re}\nabla’^2\vec{u}’-\frac{1}{Re}\frac{\epsilon}{Da}\vec{u}’-\frac{\epsilon C_F|\vec{u}’|}{\sqrt{Da}}\vec{u}'[/latex]

至于连续方程，在无量纲nabla算子下也是差个常数，仍然为：

[latex]\nabla’\cdot\vec{u}’=0[/latex]

无量纲坐标的特点是，在该坐标系下，那个特征尺度为1.不是网格雷诺数。

最后是 Courant 数。这是一个无量纲参数。其定义为[latex]C=U\Delta t/\Delta x[/latex]，也就是[latex]C=\Delta t’/\Delta x'[/latex].

也有一种姿势是[latex]C=u\Delta t/\Delta x[/latex]，u是网格中最大速度什么的，那么就是[latex]C=u’\Delta t’/\Delta x'[/latex]

参考文献

[chen2008] Chen X, Yu P, Winoto S H, et al. Numerical analysis for the flow past a porous square cylinder based on the stress-jump interfacial-conditions[J]. International Journal of Numerical Methods for Heat & Fluid Flow, 2008, 18(5): 635-655.

[nithiarasu2002] Nithiarasu P, Seetharamu K N, Sundararajan T. Finite element modelling of flow, heat and mass transfer in fluid saturated porous media[J]. Archives of Computational Methods in Engineering, 2002, 9(1): 3.