UnrealEngine 4/5 学习笔记：Niagara 流体模拟（1）
流体模拟概念流体模拟就是通过算法生成数据来呈现气体或液体等流体的运动。该模拟数据可以表示为网格或粒子，具体取决于所使用的算法。虚幻引擎流体模拟系统使用网格来模拟气体，并混合使用粒子和网格来模拟液体。网格2D集合（Grid 2D Collection） 和 网格3D集合（Grid 3D Collection） 是分别用于在 2D 和 3D 网格中存储命名属性的数据接口，这些接口用于解算流体时所需的所有计算。
模拟具有逼真移动效果的流体时，一个主要组成部分是"压力解算"过程。该技术涉及求解方程组，该方程组旨在确保流体正确地围绕对象流动，并产生逼真的漩涡运动。
气体模拟用网格表示，其中每个单元都包含表示该位置介质密度、温度和速度的数据。网格单元的尺寸越小，模拟的质量越高，但是质量提高的同时，计算成本会增加。在渲染烟雾模拟时，一般是将密度网格可视化。密度越高越不透明。火模拟的功能类似，其中温度控制每个网格单元中火的颜色。在火模拟中，由于浮力，温度越高会导致气体上升得更快。
流体模拟基础
矢量微积分
物理学中经常需要研究和解释矢量场（如电场、磁场、流场等）的物理现象问题，为了更好地解释各种“场”的性质，引入了矢量微积分的数学方法，通过对矢量场求散度、旋度可以用于判断矢量场是否有源、是否有旋，从而区分不同的场的性质和区别。如静电场是无旋的有源场，磁场是无源的涡旋场，而流场更为复杂一些，既有有旋流动又有无旋流动。
梯度（Gradient）
函数 F(x,y,z) 的梯度： grad F=\nabla F=\frac{\partial F}{\partial x}\vec{i}+\frac{\partial F}{\partial y}\vec{i}+\frac{\partial F}{\partial z}\vec{k}
向量微分（Nabla）算子： \nabla=\frac{\partial }{\partial x}\vec{i}+\frac{\partial }{\partial y}\vec{j}+\frac{\partial }{\partial z}\vec{k}
梯度是向量（矢量），表示的是某一函数在该点处的方向导数沿着该方向取得最大值，即函数在该点处沿着该方向（此梯度的方向）变化最快，变化率最大（为该梯度的模）。对于多元函数 F(x_1,\cdots,x_n)，梯度的表现方式则是向量场，这个场描述的是这个多元函数在空间中每一个点处变化的大小和方向。
散度（Divergence）
函数 F(x,y,z) 的散度： div \vec{F}=\nabla · \vec{F}=\frac{\partial F_{x}}{\partial x}+\frac{\partial F_{y}}{\partial y}+\frac{\partial F_{z}}{\partial z}
散度是标量，表示的是一个闭合曲面内单位体积的通量，用于衡量在某一点处相应的矢量聚集或者发散程度，测量方向为径向。以气体举例，空间中某一个极小的闭合区域散度就是在此区域中的气体是否有‘逸出’，逸出气体的体积就是散度的大小。散度为 0 的场，称为无源场。
旋度（Curl）
函数 F(x,y,z) 的旋度： curl \vec{F}=\nabla \times \vec{F}=(\frac{\partial F_{z}}{\partial y}-\frac{\partial F_{y}}{\partial z})\vec{i}+(\frac{\partial F_{x}}{\partial z}-\frac{\partial F_{z}}{\partial x})\vec{j}+(\frac{\partial F_{y}}{\partial x}-\frac{\partial F_{x}}{\partial y})\vec{k}
旋度是向量，表示的是单位面积的环量（环量面密度），用于衡量围绕空间中某一点的旋转速度，测量方向为切向。旋度为 0 的场，称为无旋场。
拉普拉斯算子（Laplacian）
拉普拉斯算子： \Delta F=\nabla^{2}F=\nabla·\nabla F=\frac{\partial^{2} F}{\partial^{2} x}+\frac{\partial^{2} F}{\partial^{2} y}+\frac{\partial^{2} F}{\partial^{2} z}
拉普拉斯算子是 N 维欧几里德空间中的一个二阶微分算子，定义为梯度的散度。由于标量函数的梯度往往是一种“驱动力”（或者叫“趋势”），所以针对“驱动力”求散度就可以知道空间中“源”的分布。偏微分方程 \nabla · \nabla F=0 被称为拉普拉斯方程；而如果右边为某个非 0 常数，即 \nabla · \nabla F=q ，我们称之为泊松方程。更一般地，如果梯度再乘上一个标量 a （如 \frac{1}{\rho} )，即 \nabla · (a\nabla F)=q ，我们依旧称之为泊松问题。
Naiver-Stokes 方程
纳维-斯托克斯方程（Navier-Stokes equation）是描述粘性不可压缩流体动量守恒的运动方程，简称N-S方程。
\begin{align} \frac{\partial \vec{u}}{\partial t}+\vec{u}·\nabla \vec{u}+\frac{1}{\rho}\nabla p&=\vec{g}+\nu\nabla·\nabla \vec{u} \tag{1}\\ \nabla ·\vec{u}&=0 \tag{2} \end{align}
\vec{u}=(u,v,w) 表示流体的速度矢量； \rho 表示流体的密度； p 表示压力（流体对任何物体施加的单位面积力）； \vec{g} 表示重力加速度； \nu 表示流体的运动粘度（衡量流体的黏滞性）。
动量方程
偏微分方程（1）被称为动量方程，描述的是施加在流体上的力影响流体运动的情况，是牛顿第二定理 \vec{F}=m\vec{a} 针对流体运行情况的变形。对施加在流体上和流体内部的作用力进行研究时，为了便于分析，我们将流体中无限小的一滴流体看作是微观的流体粒子，在经典物理中，每一个微观粒子都应该服从牛顿运动定律，所以流体粒子所受到的力包括重力 m\vec{g} 、压力 V\nabla p 和粘滞力 Vμ\nabla\nabla \vec{u} 。综合流体粒子所受到的力以及流体粒子加速度的定义 \vec{a}=\frac{D\vec{u}}{Dt} 可得：
\vec{F}=m\vec{g}-V\nabla p+V\mu\nabla ·\nabla \vec{u}=m\frac{D\vec{u}}{Dt} \tag{3}\\ V 表示流体粒子的体积； \mu 表示动力粘度系数； D 表示物质导数（对流体质点求导数）。
引入粒子密度 \rho=m/V 和运动粘度 \nu=\mu/\rho 对方程（3）进行简化可得：
\frac{D\vec{u}}{Dt}+\frac{1}{\rho}\nabla p=\vec{g}+\nu\nabla · \nabla \vec{u} \tag{4}\\拉格朗日描述 & 欧拉描述
研究流体的运动时，通常用两种方法来描述：拉格朗日描述（ Lagrangian viewpoint）、欧拉描述（Eulerian viewpoint）。
拉格朗日描述是指观察空间中具体的一个流体粒子，关注该流体粒子在任意时间下的运动参数，描述的是各个流体粒子的运动参数（位置坐标、运动速度、加速度等）随时间 t 变化的规律。综合所有的流体粒子的运动参数便可得到整个流体的运动规律，常用于研究波动问题。假设某一流体粒子在时间 t=0 时刻所处的位置坐标是 (a,b,c) ，则可得该流体粒子的运动方程为 \vec{x}=\vec{x}(a,b,c,t) ，流体粒子的速度 \vec{u}=\frac{\partial\vec{x}(a,b,c,t)}{\partial t} ，流体粒子的加速度为 \vec{a}=\frac{\partial \vec{u}}{\partial t}=\frac{\partial^{2}\vec{x}(a,b,c,t)}{\partial^{2}t} 。
欧拉描述是指观察空间中某一个固定点，关注该固定点上的流体物理属性，描述的是各个固定点上物理属性（如密度、速度、温度等）随时间 t 变化的规律。由于在整个观察期间，流经固定点的流体粒子会有许多个，不同流体粒子在经由固定点的前后会有加速度的变化，所以在描述固定点的速度时位置坐标 (a,b,c) 与时间是有关的。假设固定点的的运动方程为 \vec{x}=\vec{x}(a,b,c,t) ，固定点的速度 \vec{u}=\frac{d\vec{x}(a,b,c,t)}{d t} ，固定点的加速度为 \vec{a}=\frac{d\vec{u}}{dt}=\frac{\partial \vec{u}}{\partial t}+\nabla \vec{u}·\frac{d\vec{x}(a,b,c,t)}{dt} 。
定义物质导数为 <span class="ztext-math" data-eeimg="1" data-tex="\frac{D\vec{u}}{Dt}=\frac{\partial \vec{u}}{\partial t}+\nabla \vec{u}·\vec{u} \tag{5}\\">\frac{D\vec{u}}{Dt}=\frac{\partial \vec{u}}{\partial t}+\nabla \vec{u}·\vec{u} \tag{5}\\
将方程（5）代入方程（4）中即可得到动量方程（1）。
不可压缩性
偏微分方程（2）被称为不可压缩条件。为了简化计算，在计算机的流体模拟中通常把所有的流体当作是不可压缩的，即它们的体积不会发生变化。满足不可压缩条件的速度场被称为是无散度的，即在该速度场下流体体积既不膨胀也不坍缩，而是保持在一个常量。任取流体的一部分，设其体积为 V ，而其边界闭合曲面为 \partial V ，我们可以通过围绕边界曲面对流体速度 \vec{u}在曲面法线方向上的分量进行积分来衡量这块部分流体的体积变化速率： \frac{dV}{dt}=\int\int_{\partial V}\vec{u}·ndS ，对于不可压缩的流体，其体积变化速率应当为：
\int\int_{\partial V}\vec{u}·ndS=0 \tag{6}\\
由高斯散度定理（流过给定区域的边界的通量等于内部所有散度的贡献），我们可以把式（6）转换为体积分： \int\int_{\partial V}\vec{u}·ndS=\int\int\int_{V}\nabla·\vec{u}dV=0 \\
于是可以推断出 \nabla·\vec{u}=0 。
模拟不可压缩流体的关键部分就是使得流体的速度场保持无散度的状态，这也是流体内部压力的来源。为了把压力与速度场的散度联系起来，我们在动量方程（1）两边同时取散度：
\nabla·\frac{\partial \vec{u}}{\partial t}+\nabla·(\vec{u}·\nabla \vec{u})+\nabla·\frac{1}{\rho}\nabla p=\nabla·(\vec{g}+\nu\nabla·\nabla \vec{u}) \\ 将上式中第一项改变微分顺序并考虑不可压缩条件可得 \nabla·\frac{\partial \vec{u}}{\partial t}=\frac{\partial}{\partial t}\nabla·\vec{u}=0 ，所以关于压力的方程为
\nabla·\frac{1}{\rho}\nabla p=\nabla·(-\vec{u}·\nabla \vec{u}+\vec{g}+\nu\nabla·\nabla \vec{u}) \\
丢掉粘度项
在大多数的计算机流体动画模拟中，粘滞力的影响微乎其微，为秉持着方程组越简单越好的原则，我们常常丢弃粘度项。这种无粘度的流体是理想流体，而丢弃了粘度项的 Navier-Stokes 方程被称为欧拉方程：
\begin{equation*} \begin{cases}\frac{D\vec{u}}{Dt}+\frac{1}{\rho}\nabla p=\vec{g} \\ \nabla ·\vec{u}=0 \\ \end{cases} \end{equation*} \\
边界条件
在流体模拟中我们仅仅关注两种边界条件：固体墙（solid walls）、自由面（free surfaces）。
固体墙是流体与固体接触的边界，流体既不会流进固体内部也不会从固体内部流出。当固体静止不动时，流体在固体墙法线方向上的分量为 0 ： \vec{u}·n=0 ；当固体以一定的速度运动时，流体速度在法线方向上的分量与固体的移动速度在法线方向上的分量则应该保持一致： \vec{u}·n=\vec{u}_{solid}·n 。
自由面是流体与另一种流体相接触的边界，从物理学上讲，两种不同的流体之间实际上存在着与表面平均曲率成正比的压力骤变： =\lambda k 。其中 记为压力之差， \lambda 是表面张力系数。
N-S方程的分步求解
用计算机对 N-S 偏微分方程进行离散化求解时，为了程序的松耦合性以及使计算尽可能地高效、简单，在流体模拟过程中·，需要将流体方程离散化。对于不可压缩的无粘度流体方程（欧拉方程），我们将其离散化成对流项（advection）（公式 7）、体积力项（body force）（公式 8）、压力/不可压缩项（公式 9）。
\begin{align} \frac{Dq}{Dt}=0 \tag{7}\\ \frac{\partial \vec{u}}{\partial t}=\vec{g} \tag{8}\\\end{align}
\begin{equation*} \begin{cases} \frac{\partial \vec{u}}{\partial t}+\frac{1}{\rho}\nabla p=\vec{g} \\ \nabla ·\vec{u}=0 \\ \end{cases}\end{equation*} \tag{9}\\ 我们记对流项公式（7）的对流计算算法为 advect(\vec u, \Delta t, q) ，即对于给定的时间步长 \Delta t 和速度场 \vec u ，对物理量 q 进行对流。对于体积力项（8），我们采用简单的前向欧拉法即可： \vec u \leftarrow \vec u + g\Delta t，对于压力/不可压缩项（9），我们用 project(\Delta t, \vec u) 算法计算出正确的压力以确保速度场 \vec u 的无散度性质。伪代码如下：
Fluid Simulation(\vec u_n , \Delta t):\\ 1.初始化速度场\vec u_n,使得\vec u_n无散度\\ 2.对于每个时间步n = 0 , 1 , 2 , . . . \\ 3.决定一个合理的时间步长\Delta t = t_{n+1}-t_n\\ 4.对流项计算\vec u_A=advect(\vec u_n,\Delta t,\vec q) \\ 5.体积力项计算\vec u_B=\vec u_A+\Delta t\vec g \\ 6.无散度投影\vec u_{n+1}=project(\Delta t,\vec u_B)
时间步长
在流体模拟算法中，确定适当的时间步长是算法的第一步。因为计算流体模拟的最后结果是呈现在屏幕上的，所以 \Delta t 的选取与屏幕的刷新率有重要的关系。若选取的 \Delta t 有 t_{n}+\Delta t>t_{frame} ，那么必须做一个截断使 \Delta t=t_{frame}-t_n 。此外，流体模拟的三个步骤即对流项、体积力项、无散度投影项对时间步长 \Delta t 的要求不尽相同，要选择一个满足所有要求的最小时间步长能确保计算的收敛性。此外，一方面为了流体模拟的真实性，我们可能需要选取一个足够小的时间步长来复现流体的高质量细节。另一方面，有时高性能的需求又使得我们不能选取太小的时间步长去渲染一帧。假设一帧至少要进行三个时间步的模拟，那么 \Delta t 应该至少设成帧间隔时间的三分之一。
网格结构
Harlow 和 Welch 提出了一种经典的 MAC（marker and cell）网格结构，许多不可压缩流体模拟的算法都在这个网格结构上呈现出了良好的效率。MAC 网格是一种交叉排列的网格，不同类型的物理量被存储于网格的不同位置。以二维的网格为例，流体粒子的压力数据存储于网格的中心点 P_{i,j} ，而速度则沿着笛卡尔坐标被分成了两部分。水平方向的 u 成分被存储在了网格单元竖直边的中心处，例如网格单元 (i,j) 和 (i+1,j) 之间的水平速度记为 u_{i+1/2,j} 。垂直方向的 v 成分则被存储在了网格单元水平面的中心上。这样的存储方案十分有利于估算流体流进/流出某个网格单元的量。扩展到三维的情况，MAC 网格同样是交错排列的结构网格。压力数值存储在立方体网格单元的中心，三个速度分量分别被记录在立方体网格单元的三个表面的中心点上。在数值计时，这样的分配方式使得我们可以准确地采用中心差分法计算压力梯度和速度的散度，同时克服了中心差分法的一个普遍的缺点。
参考
[*]^https://docs.unrealengine.com/5.0/zh-CN/fluid-simulation-in-unreal-engine---overview/
[*]^https://yangwc.com/2019/05/01/fluidSimulation/
[*]^https://zhuanlan.zhihu.com/p/97545154
[*]^https://www.zhihu.com/question/26822364
[*]^https://blog.csdn.net/qq_31615919/article/details/86669522
[*]^https://zhuanlan.zhihu.com/p/197473887
页:
[1]