基于位置约束的方法

1. 简介

基于位置约束的方法，即PBF（Position Based Fluid）法，是一种快速、高效、非物理的不可压缩流体仿真方法。该方法通过直接对流体粒子施加恒定密度约束来更新粒子的位置与速度，进而实现不可压缩性流体运动的模拟。

2. 方法

流体在运动过程中需要时刻保持密度恒定，即恒定密度约束。粒子$i$的恒定密度约束的具体形式为：$\rho_i = \rho_0$。其中$\rho_0$为流体的静止密度，在PeriDyno中该静止密度的取值一般为$1000$；$\rho_i$为粒子$i$的密度，可通过基于核函数的SPH法插值计算获得，即：

$$ {\rho_i = \sum_{j,r_{ij}<h}m_j \cdot W(r_{ij},h)} $$ (公式1)

其中，$r_{ij}$为粒子$i$与其邻域粒子$j$的间距，即：$r_{ij} = |\mathbf{x}_i- \mathbf{x}_j|$。因此粒子$i$的密度可以被视为粒子$i$及其邻域粒子$j$空间坐标作为变量的函数。因此可获得以流体粒子空间位置作为变量的恒定密度约束条件，即： $$ C_i(\mathbf{x}) = \frac{\rho_i(\mathbf{x})}{\rho_0}-1=0 $$ (公式2)

其中，$\mathbf{x}$为包含有粒子$i$及其邻域粒子$j$空间位置坐标的矢量。

若流体中粒子$i$位置上的不可压缩性条件被破坏：$\rho_i > \rho_0$，也即是 $C_i(\mathbf{x}) > 0$ 则需要修正粒子$i$及其邻域粒子$j$的空间位置坐标，以保持粒子$i$位置附近的不可压缩性。

$$ C_i(\mathbf{x}+\triangle \mathbf{x}) = \frac{\rho_i(\mathbf{x +\triangle \mathbf{x}})}{\rho_0}-1=0 $$ (公式3)

其中$\triangle \mathbf{x}$为粒子$i$及其邻域粒子$j$的空间位置修正矢量。

粒子$i$及其邻域粒子$j$的空间位置修正矢量$\triangle\mathbf{x}$可以通过泰勒近似和牛顿法获得。略去推导过程，$\triangle\mathbf{x}$的计算式为：

$$ {\triangle \mathbf{x_{i}}}= {\frac{1}{\rho_0}}\sum_{j,r_{ij}<h}(\lambda_j+\lambda_i)\nabla W(r_{ij}, h)， $$ $$ \nabla W(r, h) = \frac{\mathrm{d} r }{\mathrm{d} \mathbf{x}} \frac{\mathrm{d} W(r ,h)}{\mathrm{d} r} $$ (公式4)

其中，$W$为SPH核函数；$\lambda$的计算式为：

$$ \lambda_i = - \frac{Ci(\mathbf{x})}{\sum_k|\nabla_k C_i|^2} $$ (公式5)

计算获得粒子空间位置修正量$\triangle \mathbf{x}$后，即可获得满足不可压缩性条件的粒子位置，按照该位置来更修正粒子位置，然后根据粒子位置修正量来更新粒子的速度，即$\mathbf{v}^{new}_i=\mathbf{v}^{old}_i + \triangle \mathbf{x_i}/{\delta t}$，来模拟流体质点的惯性作用，即可完成当前帧流体粒子运动的模拟计算。

一帧中仅计算一次粒子的位置修正量$\triangle \mathbf{x}$并不能很好地消除密度误差，为了获得更为精确的计算结果，常需要在一帧中多次迭代计算$\triangle \mathbf{x}$。

3. 实现

PeriDyno中，PBF求解器的源文件为：

"Engine/Dynamics/Module/DensityPBD.cu";
"Engine/Dynamics/Module/DensityPBD.h";

PBF求解器的输入为不满足恒定密度约束的流体粒子的位置，输出为满足恒定密度约束的粒子位置，以及惯性作用下的粒子速度。

a. 包含的模块：

  • SummationDensity（密度计算模块）

      std::shared_ptr<SummationDensity<TDataType>> m_summation;

    DensityPBD求解器使用该模块用于计算粒子的质量密度，并以此来计算恒定密度约束条件下粒子位置的修正值。

  • SpikyKernel (核函数)

     SpikyKernel<Real> m_kernel;

    该模块即为SPH法的核函数$W(r_{ij},h)$。

b. 变量：

  • 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}$)。

  • 内部变量

          DArray<Real> m_lamda;
      DArray<Coord> m_deltaPos;

    m_lamda 为上述公式5中的 $\lambda$； m_deltaPos为公式4中的 $\triangle\mathbf{x}$；

c. 参数：
    int IterationNumber = 5;
    Real RestDensity = 1000.0f;
    Real SamplingDistance = 0.005f;
    Real SmoothingLength = 0.01f;

上述参数的含义分别为：

  • IterationNumber： 迭代次数（迭代次数越多，计算结果越精确，每一帧的密度误差越小）；
  • RestDensity： 静止密度（即公式2中的$\rho_0$）；
  • SamplingDistance： 粒子的初始间距；
  • SmoothingLength： SPH法中的平滑距离/支持域半径，同时也是粒子的邻域半径。
d. 函数：
  • K_ComputeLambdas： 该函数为公式5，用于计算 $\lambda$。
  template <typename Real, typename Coord>
	__global__ void K_ComputeLambdas(...)
  • K_ComputeDisplacement： 该函数为公式4，用于计算 $\triangle \mathbf{x}$。
  	template <typename Real, typename Coord>
	__global__ void K_ComputeDisplacement(...)
  • DP_UpdateVelocity： 该函数为用于计算惯性作用下粒子的速度，即：$\mathbf{v}^{new}_i=\mathbf{v}^{old}_i + \triangle \mathbf{x_i}/{\delta t}$。
  template <typename Real, typename Coord>
	__global__ void DP_UpdateVelocity(...)

参考文献

[1] Macklin M , M Müller. Position based fluids[J]. Acm Transactions on Graphics, 2013, 32(4):1-12.