粘性是流体的重要特性。当流体产生了形变的趋势时(速度的梯度:$\nabla \mathbf{v}$),流体的粘滞力就会产生并阻滞流体发生形变。在粒子系统中,粘性还有另外一项功能,即是削弱流体个别粒子过于剧烈、高频的运动,使流体粒子具有相一致的运动,进而提高仿真稳定性。
PeriDyno中包含了一种改进与XSPH法的人工粘性求解器。虽然该粘性求解器并不是基于物理的求解器,但它能够非常稳定地够模拟大粘性流体,而且能够有效提高粘性流体的仿真稳定性,具有十分优异的计算效率。
比较经典的人工粘性求解算法是XSPH法,该方法通过削弱相邻流体粒子之间的速度差值来保证粒子的运动速度保持接近,来使流体的速度场分布趋于平滑。其计算公式如下所示。
$$ \mathbf{v_i}^{new} = \mathbf{v_i}^* + \epsilon \sum_{j, r_{ij}<h} \frac{m_j}{\rho_j} \left (\mathbf{v_j}^*-\mathbf{v_i}^*\right )W(r_{ij},h) $$
(公式1) 然而该方法存在较为严重的问题,即控制粘性强度的稀疏$\epsilon$的取值不能大于1,因此该方法只能用于提高低粘性流体的稳定性,而无法模拟大粘性流体。
因此PeriDyno在原XSPH人工粘性算法基础上进行了改进,使其能够稳定模拟大粘性流体。改进后的人工粘性计算公式如下。
$$
\mathbf{v_i} ^{new}=\frac{1}{b+1}\left ( \mathbf{v_i}^* + b\sum_{j, r_{ij}<h}\mathbf{v_i}^*W\left ( r_{ij},h \right ) \right )
$$
(公式2) 其中$b$的计算式为:
$$ b = \epsilon\frac{\delta t }{h} 。 $$
(公式3)
为了保证仿真的稳定性,需要在一帧中多次迭代计算公式2。
PeriDyno中,人工粘性求解器的源码文件为:
"Engine/Dynamics/Module/ImplicitViscosity.cu";
"Engine/Dynamics/Module/ImplicitViscosity.h";
人工粘性求解器的输入输出均为流体粒子速度。但相比于输入,输出的粒子速度在空间中的分布更为平滑。
NeighborList (邻域粒子列表)
DEF_ARRAYLIST_IN(int, NeighborIds, DeviceType::GPU, "");
为了实现粒子各项物理量及其微分量的插值计算,计算过程需要使用到粒子的邻域列表。该邻域列表需要在每一帧计算开始时预先由邻域查找模块生成,邻域查找模块的位置为:
"Engine\Framework\Topology\NeighborPointQuery.h"
"Engine\Framework\Topology\NeighborPointQuery.cu"
Velocity, Position (粒子速度与位置)
DEF_ARRAY_IN(Coord, Velocity, DeviceType::GPU, "Input particle velocity");
DEF_ARRAY_IN(Coord, Position, DeviceType::GPU, "Input particle position");
Velocity为粒子的速度($\mathbf{v}$);Position为粒子空间位置 ($\mathbf{x}$)。
DEF_VAR(Real, Viscosity, 0.05, "");
DEF_VAR(int, InterationNumber, 3, "");
DEF_VAR_IN(Real, SmoothingLength, "");
template<typename Real, typename Coord>
__global__ void IV_ApplyViscosity(...)
该函数的功能是计算公式2,获得粘性作用下的粒子速度。
Schechter H, Bridson R. Ghost SPH for animating water[J]. ACM Transactions on Graphics (TOG), 2012, 31(4): 1-8.