物理模拟——PBF流体模拟（1）

本文是有关我大学毕设的一个总结，毕设题目为：基于粒子法流体动力学的物理仿真引擎开发，实际工作为基于图形API，搭建起一套简单的渲染框架，并且基于此框架实现CPU端流体模拟算法Position Based Fluid，仓库地址为：https://gitee.com/FlyingZiming/fluid-simulation-engine

目录

• 介绍
• 基于物理的动画
  • 基础概念
  • 流体动画
• 流体模拟基础
  • 梯度
  • 散度
  • 旋度
  • 拉普拉斯算子
  • 粘性
  • 纳维-斯托克斯偏微分方程组
  • 拉格朗日视角和欧拉视角
• SPH
  • 原理及公式
  • 过程
  • 加速
• 引用

介绍

PBF全称Position Based Fluid，意为基于位置的流体，是基于Position Based Dynamics实现的流体模拟，本文会讲解物理动画概念、流体模拟基础、SPH、PBD、PBF以及流体渲染。

笔者水平有限，有说的不对的地方欢迎指出。

基于物理的动画

基础概念

基于物理的动画意味着利用物理原理，通常是经典力学，去为复杂的现象建模，模拟时间推进，实现成模拟算法。

他们的目标是使用一系列的建模方法，随着时间推演，生成一系列的系统状态，再把这些状态渲染为图像。这种动画方法通常是作用于复杂的动画，其表现通常是不被预设好的，根据环境变化有所变化的，因此无法按照传统的制作关键帧方式来进行动画制作，比如一个弹性球落地弹起、布料随风飘动、雨水落在地面。

创建物理模拟动画有两个元素：模型和模拟。模型是一系列定律或规则，通常是物理定律。模拟是将模型封装进一个框架之中从而预测该模型在时间推演下的演化。基于物理的动画的模型就是物理的，是模拟现实世界的行为，但这种模型其实也是被简化过的，影响较小的部分通常会被忽略，模型注重宏观行为，在具体细节上不是特别在意。

具体的基于物理的动画在使用中通常要考虑很多要素，动画希望动画具有实时性好、相对稳定的以及可以方便控制等等的性质。为了达到这些效果，动画师有机会重新发明物理定律，以不同的方式描述我们的世界，只要他们模拟的时候看起来像。

流体动画

3d游戏中的流体动画一般的实现比较复杂，不同规模的流体有不同的做法，甚至不同的材料都会使用不同的模拟方法。

小规模的流体，比如一杯水，可以使用SPH等粒子法去模拟，游戏中典型的实时粒子法就是Position Based Fluid，NVIDIA 的Flex便是使用此算法模拟的流体，跑在GPU上的；

中等规模的流体，通常使用网格法，欧拉视角求解NS方程，网格法适合处理时常充满空间的流体，比如烟雾等；

大规模的流体，一般就不能使用物理的方法了，而是使用波函数叠加的方式，去模拟一个场，用在海平面的高度场模拟，具体如Simulating Ocean Water。

流体模拟基础

接下来简单讲解一下针对模拟领域需要使用的基础知识、概念等，不需要完全理解。

梯度

梯度是一个矢量，他的每一个分量是一个关于空间一个轴的偏导数，2维的梯度见式：

\[\nabla f(x,y) =(\frac {\partial f}{\partial x},\frac {\partial f}{\partial y}) \]

3维的梯度见式：

\[\nabla f(x,y,z) =(\frac {\partial f}{\partial x},\frac {\partial f}{\partial y},\frac {\partial f}{\partial z}) \]

形象的理解，梯度就是函数在空间变化最快的方向。梯度可以用来估算函数的值，具体表现为同一阶泰勒泰勒展开一样的表现形式。

散度

散度面向的对象是一个向量场，用于衡量某一点的矢量的聚集或者发散程度，结果是一个标量可以将其理解为一个矢量和梯度算子的点乘，值为该矢量的三个空间偏导数的和，如果矢量场的散度为0，称该矢量场无散度。

二阶的散度表示见式：

\[\nabla \cdot \vec{u} = \nabla \cdot (u,v) = \frac{\partial u}{\partial x}+ \frac {\partial v}{\partial y} \]

旋度

旋度表示围绕某一点的旋转速度，三维形式旋度是一个向量，二维形式的旋度是一个标量，旋度不能简单的进行维度推广，旋度符号可以理解为是梯度算子和一个矢量的叉乘。二维的旋度表示见式：

\[\nabla \times \vec{u} = \nabla \times (u,v) = \frac{\partial v}{\partial x}- \frac {\partial u}{\partial y} \]

拉普拉斯算子

拉普拉斯算子定义为梯度的散度，符号表示为 ，拉普拉斯算子是一个二阶的微分算子，偏微分方程 被称为拉普拉斯方程，右边如果是非零项则被称为泊松方程。拉普拉算子2维的表示见式：

\[\nabla \cdot \nabla f =\frac{\partial^2 f}{\partial x^2}+ \frac {\partial^2 f}{\partial y^2} \]

粘性

粘性是流体持续剪切变形时内部产生剪切力的性质，他宏观上表现为流体之间的动摩擦力，本质上是分子间的作用力，不同速度分子的化学键的断裂和重新生成，两层液体分子之间的牵扯和挤压。

描述粘性使用两种粘性系数，动力粘性系数和运动粘性系数。

$\mu $ 动力粘性系数表示相同变形率是的粘性大小， $ \upsilon $ 运动粘性系数表示的是相同加速度时的粘性大小，在Navier-Stokes方程中出现的即是运动粘性系数。

纳维-斯托克斯偏微分方程组

流体模拟面临的最重要的方程就是Navier-Stokes方程，不论以何种模拟方法，都是求解这个流体运动方程，见式：

\[\frac{\partial \vec{u}}{\partial {t}} + \vec{u}\cdot\nabla\vec{u}+\frac{1}{\rho}\nabla p = \vec{g}+\upsilon \nabla \cdot \nabla \vec{u} \]

\[\nabla \cdot \nabla \vec{u} = 0 \]

该方程表示的是流体的动量方程，本质上是一个牛顿第二定律的形式，描述流体上的力是如何影响流体的运动的，流体粒子所受的力包括粒子的重力、粒子之间所给与的压力、流体粒子之间相互运动产生的粘滞阻力。

Navier-Stokes方程是一个二阶非线性偏微分方程，求解起来非常复杂，在模拟领域使用计算机对NS方程进行求解，通常将方程分解为几个部分，然后按顺序推进。

拉格朗日视角和欧拉视角

拉格朗日视角和欧拉视角是研究流体和可变形固体运动的两大常用方法，其中拉格朗日视角一般对应粒子法，而欧拉视角一般对应网格法，近年来也出现了混合欧拉拉格朗日视角的方法。

拉格朗日视角基于粒子，将流体的数据，比如密度、质量、速度全部都记录在一个具体的粒子上，粒子会随着模拟的推进而改变位置，这个粒子就相当于一个数据采样点，更加具体一点理解，可以将每个粒子当成一个流体微团，拉格朗日方法模拟的是流体微团之间的相互作用得到的运动。这是一种对物质本身进行离散的NS方程求解方法，这种方法有很多的好处，比如：理解起来直观、模拟中的Advection部分是易于守恒的，但是拉格朗日视角下有不擅长的点，其中最大的就是粒子的邻居查找，这块需要使用空间的数据结构，每次执行算法前都需要更新粒子的邻居结构，这块查找邻居的时间一般是算法的瓶颈。

欧拉视角基于网格，将流体的数据记录在一组固定的空间点上面，这个空间点上的属性的含义表示通过这个空间点的属性是多少，空间点位置是不会随着模拟的推进而改变的。这是一种对空间进行离散的NS方程求解方法，这种做法的好处就是他在处理整个场的数据是自然离散到了空间点上，在做邻居查找时也十分的遍历，可以通过固定的点加上偏移得到周围空间点的数据，但是欧拉视角在处理粒子运动时经常会出现能量耗散现象，导致整个流体看起来粘性很大，他并不想拉格朗日视角处理Advection那样的自然的守恒。

两视角的主要区别就在于他们的采样点是否随着模拟推进而运动，胡渊鸣说：拉格朗日视角是随波逐流，欧拉视角是岿然不动。

近年来最火的是MPM物质点法，他属于混合欧拉拉格朗日方法，但是笔者没有具体了解过所以就不多说了。

SPH

原理及公式

Smoothed-Particle Hydrodynamics光滑粒子流体动力学是一种用来模拟连续介质动力学的离散方法，他属于拉格朗日法，最初用于天体物理领域，近年来在计算机图形学领域热度逐渐提高，常用来模拟流体运动。由于其无网格的特性，SPH在模拟复杂边界、扭曲拉伸、湍流上具有天然的优势。

SPH的高层思想就是使用一堆随着模拟运动而更新位置的粒子去携带物理量，使用核函数来近似一个连续的场，某一个位置的物理量的值等于周围粒子的物理量的加权平均，核函数就是用来决定权值，距离该位置越近的权值越大。

变种非常多，比如：PBF、PCI-SPH、WCSPH等等，这里以WCSPH的算法作引入：

其核心的表述求一点的属性的公式为，A表示一个抽象属性，可以是密度、压力、速度等物理量：

\[A(x) = \sum _i A_i \frac{m_i}{\rho_i}W(||x-x_j||_2,h) \]

除了属性的值还有计算属性的梯度的方法，这一点上他的梯度不是由上面的公式推导而来的，而是自己做的一个近似公式，但这条公式是对称的，使用时可以保证动量守恒，所以大家还是接受这条求梯度的公式：

\[\nabla A_i =\rho_i \sum _j {(\frac{A_i}{\rho_i^2}+\frac{A_j}{\rho_j^2}) \nabla_{x_i} W(||x_i-x_j||_2,h))} \]

有了这两条公式还需要流体相关的公式，需要计算流体粒子的加速度，还有模拟受到的各种粘性力、表面张力等物理定律描述的力，粒子的加速度使用材料导数，材料导数包括了两项，粒子速度随时间的变化和粒子速度在空间上的变化，在流体力学领域速度的材料导数表示为：

\[\frac{Dv}{Dt} = -\frac{1}{\rho}\nabla p +g \]

加速度公式涉及到压强梯度、密度以及外力g。

计算密度就是用上面的抽象属性公式，将密度代入A：

\[\rho(x) = \sum _i {m_i}W(||x-x_j||_2,h) \]

计算压强的公式使用了Tait方程，是一个理想气体状态方程的变种，不过实际使用还是用的SPH梯度算法来求的

\[p = B((\frac{\rho}{\rho_0})^\gamma - 1) \]

B：bulk modulus，$ \rho_0$ 是一个理想密度，密度很大压强就很大，压强很高就会尝试推着周围的粒子离开自己。

过程

有了这些公式，现在就可按照一个模拟的积分公式来确定整个SPH的算法流程：

1、对每个粒子计算它所在位置的密度；

2、对每个粒子计算它的压强梯度，进而得到它的Dv/Dt（其实就相当于求一个力）；

3、Symplectic Euler Step半隐式欧拉方法推进模拟，先根据算出来的力更新速度，然后用更新过的速度来更新位移。

如此一来SPH的流程遍完成了，这是WCSPH的做法，半隐式的时间积分经常需要考虑CFL条件问题，从材料刚度和粒子运动速度方面考虑有时间步长的限制，CFL条件是一个简单而直观的判断计算是否收敛的必要条件。它的比较直观的物理解释就是时间推进求解的速度必须大于物理扰动传播的速度，因为只有这样才能将物理上所有的扰动俘获到。不满足的话整个模拟容易发生数值爆炸等问题，满足CFL条件意味着当单位时间位移值和时间步长趋于取极限时，数值计算所求的解就会收敛到原微分方程的解。（有关CFL条件，笔者现在也不是理解很透彻，在代码里，他主要是使得我们的代码是变步长时间积分的，而不是定步长的）

加速

SPH算法在使用过程中一般最大瓶颈就在粒子的邻居搜索，因为粒子随模拟运动而运动，邻居在每次算法执行前都需要进行更新，要精确的找到一个粒子的周围粒子这个直接暴力的时间复杂度是O（n^2），为了降低这个瓶颈带来的影响加速SPH算法，一般会使用特殊的数据结构来存储粒子邻居，比如一个将位置信息通过哈希算法投射到一个数组中，一般的步骤是先将连续的粒子位置信息离散为几个大块，分成格子，然后做格子位置到粒子数组的映射，格子之间访问邻居非常简单只需要加上偏移量。这样一来一个粒子周围的粒子就等于他所在的格子的所有粒子加上周围一圈格子的所有粒子，每次算法结束步骤时进行维护，根据位置再对格子里的粒子做一个更新。

SPH邻居搜索的加速方法是很多，比如使用CPU并行进行邻居搜索、使用CUDA等GPU计算来加速。

引用

https://yangwc.com/2019/05/01/fluidSimulation/

https://diglib.eg.org/bitstream/handle/10.2312/PE.PG.PG2012short.029-034/029-034.pdf

https://book.douban.com/subject/35287308/

2021.8.25