【深入浅出 Nvidia FleX】(2) Position Based Fluids

【深入浅出 Nvidia FleX】(2) Position Based Fluids

上一篇文章 【深入浅出 Nvidia FleX】(1) Position Based Dynamics 介绍了Nvidia FleX的理论基础PBD方法。PBD一开始是为模拟固体而设计的,Miles等人[1]将它与SPH插值方法结合,提出了基于PBD框架的流体模拟方法Position Based Fluids(PBF),这也是FleX中流体模拟方法的理论基础。


一、PBF

  1. 不可压缩性约束与拉格朗日乘子推导

在文章 电影工业中的流体模拟(七)----PCISPH 中,我们介绍了PCISPH方法。 其中,流体粒子 i 的密度由SPH插值计算而得:

\rho_i = \sum\limits_j m_j W({p}_i-{p}_j, h) (1)

这里, m_j 是邻居粒子 j 的质量。 p_ip_j 为流体粒子 ij 的位置, h 为光滑核半径。

在不可压缩流体模拟中,我们希望粒子 i 的密度尽量与静态密度 \rho_0 相同,即 \rho_i = \rho_0 。因此,需要针对每一个流体粒子都施加一个常量密度约束,PBF将该约束定义为:

C_i({p}_1,...,{p}_n) = \frac{\rho_i}{\rho_0} -1 (2)

将公式(1)代入(2)中,可得:

C_i({p}_1,...,{p}_n) = \frac{\sum\limits_j m_j W({p}_i-{p}_j, h)}{\rho_0} -1 (3)

PBF的密度计算公式(1)用了[2]的Poly6核函数:

\begin{equation} W_\text{poly6}(\mathbf{r},h) = \frac{315}{64\pi h^9} \begin{cases} (h^2 - |\mathbf{r}|^2)^3 & 0 \le |\mathbf{r}| \le h \\ 0 & \text{otherwise} \end{cases} \end{equation} (4)

而梯度 \nabla W({p}_i-{p}_j, h) 则用了[2]的Spiky核函数

\begin{equation} W_\text{spiky}(\mathbf{r},h) = \frac{15}{\pi h^6} \begin{cases} (h-|\mathbf{r}|)^3  & 0 \le |\mathbf{r}| \le h \\ 0 & \text{otherwise}, \end{cases} \end{equation} (5)

\begin{equation} \nabla W_\text{spiky}(\mathbf{r},h) = -\frac{45}{\pi h^6} \begin{cases} (h-|\mathbf{r}|)^2 \frac{\mathbf{r}}{|\mathbf{r}|}  & 0 \le |\mathbf{r}| \le h \\ 0 & \text{otherwise}, \end{cases} \end{equation} (6)

因此,粒子 i 的约束函数(3)是一个关于其自身位置 \mathbf{p}_i 及其邻居粒子 j 位置 \mathbf{p}_j 的非线性方程 C_i({p}_1,...,{p}_n)  = 0 ,所有粒子 i 的约束组成了一个非线性方程组。

:PBF只考虑粒子质量相同的情况,故在下面的推导过程中可以省去质量 m_j ,即:

\rho_i = \sum\limits_j W({p}_i-{p}_j, h)

则约束函数(3)可以写为:

C_i = \frac{\sum\limits_j W({p}_i-{p}_j, h)}{\rho_0} -1 (7)

其梯度 \nabla C_i 为:

\nabla_{{p}_k}C_i =  \frac{1}{\rho_0}\sum\limits_j\nabla_{{p}_k} W({p}_i-{p}_j, h) (8)

这里有两种情况:粒子 k 为粒子 i 自身( k = i )或粒子 i 的邻居粒子 jk = j ),则有:

\begin{equation} \nabla_{{p}_k}C_i = \frac{1}{\rho_0} \begin{cases} 	\sum\limits_j \nabla_{{p}_k} W({p}_i-{p}_j, h) & \text{if }  k = i \\ 	-\nabla_{{p}_k} W({p}_i-{p}_j, h)  & \text{if } k = j \end{cases} \end{equation} (9)

可以这么理解:

  1. k = i 时,\nabla_{{p}_k}C_i = \nabla_{{p}_i}C_i 表示约束函数 C_i 关于 {p}_i 的梯度,方向为 \frac{ \mathbf{p}_i - \mathbf{p}_j }{| \mathbf{p}_i - \mathbf{p}_j |}
  2. k = j 时, \nabla_{{p}_k}C_i = \nabla_{{p}_j}C_i 表示约束函数 C_i 关于 {p}_j 的梯度,方向为-\frac{ \mathbf{p}_i - \mathbf{p}_j }{| \mathbf{p}_i - \mathbf{p}_j |}

上一篇我们介绍过PBD方法中需要求拉格朗日乘子:

\begin{equation} \lambda =  - \frac{C(\mathbf{p})}{|\nabla_\mathbf{p} C(\mathbf{p})|^2}   \end{equation} (10)

将公式(9)代入公式(10)可得:

\begin{equation}  \lambda_i = -\frac{C_i({p}_1,...,{p}_n)}{\sum_k\left|\nabla_{{p}_k}C_i\right|^2} \end{equation} (11)

这个拉格朗日乘子  \lambda_i 对于一个约束 C_i({p}_1,...,{p}_n) 中的所有粒子而言都是一样的。


2. 软约束(Soft constraint)与混合约束力(Constraint Force Mixing, CFM)

如果一个约束条件不能被违反,我们称之为硬约束(hard constraint);反之,能够一定程度上被违反的约束为软约束(soft constraint)。理想情况下,我们当然希望约束始终都是硬约束,但是由于计算误差或者数值稳定性等原因,我们有时也需要约束呈现出soft性质。

具体到PBF,当 |\mathbf{r}| = h ,即粒子 j 处于粒子 i 的光滑核半径 h 边缘时(此时粒子 i 和粒子 j 处于分离状态),公式(6) Spiky核函数的梯度值 \nabla W_\text{spiky}(\mathbf{r},h)  = 0 。如果粒子之间都处于这种状态,则由公式(9)可知:

\nabla_{{p}_j}C_i = -\frac{1}{\rho_0}\nabla_{{p}_j} W(\mathbf{r}, h)   = 0

\nabla_{{p}_i}C_i = \frac{1}{\rho_0}\sum_j{\nabla_{{p}_i}W(\mathbf{r}, h) } = 0

从而公式(11)的分母 {\sum_k\left|\nabla_{{p}_k}C_i\right|^2} = 0 ,这会导致除零错误。为了解决这个问题,PBF借鉴了Open Dynamics Engine中的混合约束力方法[2],使密度约束(2)变成软约束。具体做法是将根据约束函数求出的约束力(PBF中拉格朗日乘子 \lambda 起到类似作用)再加入到原始约束函数中:

\begin{equation} C({p} + \Delta{p}) \approx C({p}) + \nabla C^T \nabla C\lambda + \epsilon\lambda = 0 \end{equation} (12)

这里 \epsilon 是松弛参数,可以由用户指定。

引入公式(12)之后,拉格朗日乘子计算公式(11)变为:

\begin{equation}  \lambda_i = -\frac{C_i({p}_1,...,{p}_n)}{\sum_k\left|\nabla_{{p}_k}C_i\right|^2 + \epsilon} \end{equation} (13)

从而可得粒子 i 在经过上述约束投影后对应的位移向量(包括自身密度约束以及邻居粒子密度约束共同作用的效果): \begin{equation}  \end{equation}

\begin{equation} \Delta{p}_{i} = \lambda_i \nabla_{{p}_i}C_i + \sum_j \lambda_j \nabla_{{p}_j}C_i \end{equation}

\begin{align} = \frac{1}{\rho_0}\sum\limits_j{\lambda_i\nabla_{p_i} W(\mathbf{r}, h) }  + (- \frac{1}{\rho_0}\sum\limits_j{\lambda_j\nabla_{p_j} W(\mathbf{r}, h) } )   \end{align}

\begin{align} = \frac{1}{\rho_0}\sum\limits_j{\lambda_i\nabla_{p_i} W(\mathbf{r}, h) }  + \frac{1}{\rho_0}\sum\limits_j{\lambda_j\nabla_{p_i} W(\mathbf{r}, h) }    \end{align}

\begin{align} = \frac{1}{\rho_0}\sum\limits_j{(\lambda_i + \lambda_j)\nabla_{p_i} W(\mathbf{r}, h) }  \end{align} (14)


3. Tensile Instability

对于采用SPH插值技术计算密度的流体模拟方法,通常需要30∼40个邻居粒子才能使密度求值结果趋于静态密度。在邻居粒子不足的情况下,会导致通过公式(1)求出的流体密度低于静态密度,由此造成压强为负数,原本粒子间的压力变为吸引力,使粒子产生不符合实际情况的凝聚,此即SPH的Tensile Instability问题在流体模拟中的具体体现,这导致的结果是流体表面的模拟给人感觉不真实。如下图所示:

[3]采用了一种人工排斥力计算模型,当流体粒子距离过近时该排斥力会使它们分开,从而避免粒子凝聚现象。当流体粒子的压强变为负数时,用该排斥力代替压力可以有效消除SPH方法的Tensile Instability问题,防止负压强导致的粒子间非自然吸引。PBF也采用了类似的方法,在公式(14)的基础上,加入排斥项(repulsive term) s_{corr}

\begin{equation} \Delta{p}_{i} = \frac{1}{\rho_0}\sum\limits_j{\left(\lambda_i + \lambda_j + s_{corr}\right)\nabla W({p}_i-{p}_j, h)} \end{equation}

其中,

\begin{equation}  s_{corr} = -k\left(\frac{W({p}_i-{p}_j, h)}{W(\Delta{q}, h)} \right)^n \end{equation}

这里, \Delta{\mathbf{q}} 表示距离粒子 i 固定距离的点,通常取 \left|\Delta{{q}}\right| = 0.1h\cdots0.3h 。另外, k 可以看做表面张力参数(因为 s_{corr} 可以产生类似于表面张力的效果,使流体表面粒子均匀分布),这里取值 k = 0.1n = 4

如下图所示:

:FleX中并没有采用这种方法,而是采用了更好的surface tension计算模型,后面会进行介绍。

4. Vorticity Confinement

PBD方法通常会引入额外的阻尼,导致整个系统的能量损耗,由此会导致本来该有的一些涡流(vortices)快速消失。与[4]类似,PBF通过vorticity confinement向系统重新注入能量:

\begin{equation} {f}_i^{vorticity} = \epsilon\left({N}\times{\omega}_i\right) \end{equation}

其中, {N} = \frac{\eta}{|\eta|}\eta = \nabla{|\omega|}_i ,而粒子 i 的旋度{\omega}_i 为:

\begin{equation} {\omega}_i = \nabla\times{v}= \sum\limits_j ({v}_j-{v}_i) \times\nabla_{{p}_j} W({p}_i-{p}_j, h) \end{equation}

Vorticity Confinement的基本思路是:通过添加体积力(body force) {f}_i^{vorticity} 的方式,在旋度粒子(可直观理解为比周围粒子旋转快的粒子, {\omega}_i 指向粒子 i 的旋转轴)处加速粒子的旋转运动,通过这种方式来保持系统的旋度。  \epsilon 用来控制Vorticity Confinement的强度。


5. XSPH Artificial Viscosity

SPH流体模拟方法中,人工粘性(Artificial Viscosity)除了可以增加模拟的数值稳定性,还可以消除非物理振荡(nonphysical oscillations)。拉格朗日流体模拟方法中,人工粘性本质上会对流体粒子的相对运动产生阻尼作用,使流体的动能转化为热能。与[5]一样,PBF使用了XSPH直接更新速度,通过这种方式来产生阻尼。

{\begin{equation} {v}^{new}_i = {v}_i + c\sum\limits_j {v}_{ij}\cdot W({p}_i-{p}_j, h) \end{equation}}

这里 c = 0.01


二、算法和实现

  1. 算法

PBF算法与上一篇介绍的经典PBD算法类似,只是经典PBD算法采用了顺序高斯-赛德尔(Sequential Gauss-Seidel,SGS)迭代求解,而SGS不容易被GPU并行化,因此基于CUDA实现的PBF Solver使用了雅克比(Jacobi)迭代方法并行求解。

具体算法如下图所示:


2. 实现

由于FleX没有开源,我们结合下面开源PBD项目中的代码进行介绍。(注:由于字数限制,删掉了部分代码引用,具体请参考源码)

InteractiveComputerGraphics/PositionBasedDynamicsgithub.com图标

该项目中FluidDemo实现了PBF中的算法:

首先,根据外力(这里只有重力)进行数值积分:

// Time integration
for (unsigned int i = 0; i < pd.size(); i++)
{ 
	model.getDeltaX(i).setZero();
	pd.getLastPosition(i) = pd.getOldPosition(i);
	pd.getOldPosition(i) = pd.getPosition(i);
	TimeIntegration::semiImplicitEuler(h, pd.getMass(i), pd.getPosition(i), pd.getVelocity(i), pd.getAcceleration(i));
}

其次,进行邻居粒子查找操作:

// Perform neighborhood search
START_TIMING("neighborhood search");
model.getNeighborhoodSearch()->neighborhoodSearch(&model.getParticles().getPosition(0), model.numBoundaryParticles(), &model.getBoundaryX(0));
STOP_TIMING_AVG;

然后,对流体粒子的密度约束进行投影操作

// Solve density constraint
START_TIMING("constraint projection");
constraintProjection(model);
STOP_TIMING_AVG;

在constraintProjection中设置了最大迭代次数

const unsigned int maxIter = 5; 

对应算法中的solverIterations。每次迭代又分三个步骤:

  • 步骤一:计算公式(11) 的拉格朗日乘子\lambda_i ,对应代码为:
#pragma omp parallel default(shared)
{
	#pragma omp for schedule(static)  
	for (int i = 0; i < (int)nParticles; i++)
	{
		Real density_err;
		PositionBasedFluids::computePBFDensity(i, nParticles, &pd.getPosition(0), &pd.getMass(0), &model.getBoundaryX(0), &model.getBoundaryPsi(0), numNeighbors[i], neighbors[i], model.getDensity0(), true, density_err, model.getDensity(i));
		PositionBasedFluids::computePBFLagrangeMultiplier(i, nParticles, &pd.getPosition(0), &pd.getMass(0), &model.getBoundaryX(0), &model.getBoundaryPsi(0), model.getDensity(i), numNeighbors[i], neighbors[i], model.getDensity0(), true, model.getLambda(i));
	}
}

:与PBF不同,FleX[7]中的密度约束被设置为单边(unilateral)约束:

C_i({p}_1,...,{p}_n) = \frac{\rho_i}{\rho_0} -1 \le 0

PBD方法对于约束的处理方式如下:

  1. 等式约束:总是进行投影操作。
  2. 不等式约束 C_i({p}_1,...,{p}_n) \ge 0 :只有在不等式约束条件不满足即 C_i({p}_1,...,{p}_n) \lt 0 时才进行约束投影操作。

所以,只有在上面的单边约束条件不满足,即 \frac{\rho_i}{\rho_0} -1 \gt 0\rho_i  \gt \rho_0 时,才进行约束投影操作。直观理解就是只有在粒子靠的比较近(流体压缩了)的情况下,才需要进行操作让粒子分开(保持流体不可压缩),如下图所示:

Density Constraint

而当粒子分开时,即  \frac{\rho_i}{\rho_0} -1 \le 0\rho_i  \le \rho_0,不等式约束条件满足,此时不需要进行约束投影,因此也就避免了Tensile Instability一节介绍的表面粒子凝聚问题。

对应的代码为computePBFLagrangeMultiplier()中的:

	// Evaluate constraint function
	const Real C = std::max(density / density0 - static_cast<Real>(1.0), static_cast<Real>(0.0));			// clamp to prevent particle clumping at surface


另外,在计算流体粒子密度的时候computePBFDensity()考虑了流体粒子 i 以及固定边界粒子 k 的贡献[6]:

\begin{equation} \label{eq:WeightedMassDensity}     \rho_i = \sum_j m_j W(\mathbf{x}_{ij}, h) + \sum_k \psi_{b_k} W(\mathbf{x}_{ik}, h) \end{equation}

其中,

\begin{equation}  \psi_{b_k} =  \rho_0 \frac{m_{b_k}}{\sum_l m_{b_l}W(\mathbf{x}_{kl}, h)} \end{equation} (15)

这里, k 表示固体粒子,质量为 m_{b_k} ,其邻居固体粒子 l的 质量为 m_{b_l}

尽管这种方法能够更准确地计算流体密度,但是由于增加了很多固定边界粒子,由此带来了额外的计算和存储开销。所以,FleX中并没有采用这种固定边界粒子的方法。但是,上式还考虑了动态固体表面粒子对流体粒子的影响,而FleX中刚体软体布料等也是由粒子组成的,我们可以采用类似的简化方案[7]:

\begin{equation} \label{eq:WeightedMassDensity}     \rho_i = \sum_{fluid} m_j W(\mathbf{x}_{ij}, h) + solidPressure * \sum_{solid} m_{k} W({x_{ik}}, h) \end{equation}

通过下面参数控制固体表面粒子对流体粒子的影响。

float solidPressure; //!< Add pressure from solid surfaces to particles

感兴趣的同学可以去FleX Github页面下载试试,打开Fluid Cloth Coupling Water场景调整solidPressure参数范围0~1可以看到如下效果:

solidPressurehttps://www.zhihu.com/video/1046150401157599232
  • 步骤二:计算公式(14)中的位移向量 \Delta{p}_{i} ,对应代码如下:
#pragma omp parallel default(shared)
{
	#pragma omp for schedule(static)  
	for (int i = 0; i < (int)nParticles; i++)
	{
		Vector3r corr;
		PositionBasedFluids::solveDensityConstraint(i, nParticles, &pd.getPosition(0), &pd.getMass(0), &model.getBoundaryX(0), &model.getBoundaryPsi(0), numNeighbors[i], neighbors[i], model.getDensity0(), true, &model.getLambda(0), corr);
		model.getDeltaX(i) = corr;
	}
}

bool PositionBasedFluids::solveDensityConstraint(
	const unsigned int particleIndex,
	const unsigned int numberOfParticles,
	const Vector3r x[],	
	const Real mass[],
	const Vector3r boundaryX[],
	const Real boundaryPsi[],
	const unsigned int numNeighbors,
	const unsigned int neighbors[],
	const Real density0,
	const bool boundaryHandling,
	const Real lambda[],
	Vector3r &corr)
{
	// Compute position correction
	corr.setZero();
	for (unsigned int j = 0; j < numNeighbors; j++)
	{
		const unsigned int neighborIndex = neighbors[j];
		if (neighborIndex < numberOfParticles)		// Test if fluid particle
		{
			const Vector3r gradC_j = -mass[neighborIndex] / density0 * CubicKernel::gradW(x[particleIndex] - x[neighborIndex]);
			corr -= (lambda[particleIndex] + lambda[neighborIndex]) * gradC_j;
		}
		else if (boundaryHandling)
		{
			// Boundary: Akinci2012
			const Vector3r gradC_j = -boundaryPsi[neighborIndex - numberOfParticles] / density0 * CubicKernel::gradW(x[particleIndex] - boundaryX[neighborIndex - numberOfParticles]);
			corr -= (lambda[particleIndex]) * gradC_j;
		}
	}

	return true;
}
  • 步骤三:更新粒子位置 {x}_i^* \Leftarrow {x}_i^* + \Delta{p}_i ,对应代码如下:
#pragma omp parallel default(shared)
{
	#pragma omp for schedule(static)  
	for (int i = 0; i < (int)nParticles; i++)
	{
		pd.getPosition(i) += model.getDeltaX(i);
	}
}

迭代完成之后,需要更新粒子速度,对应代码:

	// Update velocities	
	for (unsigned int i = 0; i < pd.size(); i++)
	{
		if (m_velocityUpdateMethod == 0)
			TimeIntegration::velocityUpdateFirstOrder(h, pd.getMass(i), pd.getPosition(i), pd.getOldPosition(i), pd.getVelocity(i));
		else
			TimeIntegration::velocityUpdateSecondOrder(h, pd.getMass(i), pd.getPosition(i), pd.getOldPosition(i), pd.getLastPosition(i), pd.getVelocity(i));
	}

然后可以进行我们之前介绍过的Vorticity Confinement操作,PositionBasedDynamics开源项目没有实现。FleX里实现了Vorticity Confinement,通过下面参数进行控制:

float vorticityConfinement; //!< Increases vorticity by applying rotational forces to particles

打开FleX DamBreak 5cm场景调整vorticityConfinement参数范围可以得到下面这样的效果:

Vorticity Confinementhttps://www.zhihu.com/video/1046509209734242304

之所以出现上面这样的效果,是因为vorticityConfinement向系统中重新注入能量,当取值过大时会让本来达到平衡状态的流体又开始运动起来。

另外,还可以应用XSPH Artificial Viscosity,对应代码:

/** Compute viscosity accelerations.
*/
void TimeStepFluidModel::computeXSPHViscosity(FluidModel &model)
{
	ParticleData &pd = model.getParticles();
	const unsigned int numParticles = pd.size();	

	unsigned int **neighbors = model.getNeighborhoodSearch()->getNeighbors();
	unsigned int *numNeighbors = model.getNeighborhoodSearch()->getNumNeighbors();

	const Real viscosity = model.getViscosity();
	const Real h = TimeManager::getCurrent()->getTimeStepSize();

	// Compute viscosity forces (XSPH)
	#pragma omp parallel default(shared)
	{
		#pragma omp for schedule(static)  
		for (int i = 0; i < (int)numParticles; i++)
		{
			const Vector3r &xi = pd.getPosition(i);
			Vector3r &vi = pd.getVelocity(i);
			const Real density_i = model.getDensity(i);
			for (unsigned int j = 0; j < numNeighbors[i]; j++)
			{
				const unsigned int neighborIndex = neighbors[i][j];
				if (neighborIndex < numParticles)		// Test if fluid particle
				{
					// Viscosity
					const Vector3r &xj = pd.getPosition(neighborIndex);
					const Vector3r &vj = pd.getVelocity(neighborIndex);
					const Real density_j = model.getDensity(neighborIndex);
					vi -= viscosity * (pd.getMass(neighborIndex) / density_j) * (vi - vj) * CubicKernel::W(xi - xj);

				}
// 				else 
// 				{
// 					const Vector3r &xj = model.getBoundaryX(neighborIndex - numParticles);
// 					vi -= viscosity * (model.getBoundaryPsi(neighborIndex - numParticles) / density_i) * (vi)* CubicKernel::W(xi - xj);
// 				}
			}
		}
	}
}

FleX中通过下面参数控制XSPH 粘性:

float viscosity;					//!< Smoothes particle velocities using XSPH viscosity

具体效果如下:

XSPH viscosityhttps://www.zhihu.com/video/1046511045098672128

再来看一个高粘度流体的粒子:

High Viscosityhttps://www.zhihu.com/video/1047098342210899968


之前我们介绍的流体内部的压力等都是在宏观尺度上(可以理解为流体粒子级别)进行度量的,液体(比如水)的模拟还存在着下面三种微观尺度的力(流体分子级别):

1. 表面张力(Surface Tension)。液体与气体相接触时,会形成一个表面层,在这个表面层内存在着的相互吸引力就是表面张力,它能使液体表面自动收缩。如下图所示:

:前面我们介绍Tensile Instability时讲过PBF采用的 s_{corr} 可以产生类似于表面张力的效果,下面将要介绍的方法是用更准确的物理模型来计算液体的表面张力。

2. 内聚力(Cohesion),是液体内部相邻分子之间的相互吸引力。表面张力就是由于液体分子间的内聚力而产生的。从上图可以看到,在液体内部,分子受到来自各个方向的分子作用相互抵消,合力为零;而自由表面(和空气接触的部分)的分子只会受到来自其他自由表面分子以及流体内部分子的作用力。

3. 附着力 (Adhesion),是指液体分子和固体分子之间的相互吸引力。如下图,如果没有附着力,水滴不会依附在松针表面。

我们先介绍一下[8]中的这三种力的模型:

  1. 内聚力 \mathbf{F}^{cohesion}_{i \leftarrow j}

\mathbf{F}^{cohesion}_{i \leftarrow j} = -\gamma m_i m_j C(|\mathbf{x}_i - \mathbf{x}_j|) \frac{\mathbf{x}_i - \mathbf{x}_j}{|\mathbf{x}_i - \mathbf{x}_j|} (16)

其中,流体粒子 i 和粒子 j 的质量和位置分别为 m_i\mathbf{x}_im_j , \mathbf{x}_j\gamma 为表面张力系数, C 为:

\begin{equation} C(r) = \frac{32}{\pi h^9}  \left\{ \begin{array}{1r} (h-r)^3 r^3 &  \text{if }  2r > h \land  r \le h \\ 	 2(h-r)^3 r^3-\frac{h^6}{64} &  \text{if} r > 0 \land 2r \le h \\ 0  & \text{if }  otherwise  \end{array}  \right. \end{equation}

2. 表面张力 \mathbf{F}^{surfaceTension}_{i \leftarrow j}

除了内聚力外,[8]还定义了一种curvature based force \mathbf{F}^{curvature}_{i \leftarrow j} , 使得:\mathbf{F}^{surfaceTension}_{i \leftarrow j}  =  \frac{2\rho_0}{\rho_i + \rho_j} (\mathbf{F}^{cohesion}_{i \leftarrow j} + \mathbf{F}^{curvature}_{i \leftarrow j} )

\mathbf{F}^{curvature}_{i \leftarrow j}  = -\gamma m_i (\mathbf{n}_i - \mathbf{n}_j)

其中,\mathbf{n}_i = \nabla c_s h 为缩放因子, c_s 为smoothed color field[9]

c_s = h\sum\frac{m_j}{\rho_j} W(|\mathbf{x}_i - \mathbf{x}_j|)

3. 附着力 \mathbf{F}^{adhesion}_{i \leftarrow k} :

\mathbf{F}^{adhesion}_{i \leftarrow k} = -\beta m_i \psi_{b_k} A(|\mathbf{x}_i - \mathbf{x}_k|) \frac{\mathbf{x}_i - \mathbf{x}_k}{|\mathbf{x}_i - \mathbf{x}_k|} (17)

这里 k 为固体粒子, \psi_{b_k} 依据公式(15)进行计算。 \beta 为附着力系数, A 为:

\begin{equation}  A(r) = \frac{0.007}{h^{3.25}}  \begin{cases} 	 \sqrt[4]{-\frac{4r^2}{h} + 6r -2h} & \text{if }  2r > h \land r \le h \\ 	 0 & \text{if } otherwise  \end{cases}  \end{equation}


PositionBasedDynamics开源项目没有实现上述三种力,下面结合UnifiedParticleFrameworkCUDA 中的代码作介绍:

Gfans/UnifiedParticleFrameworkCUDAgithub.com图标
void UnifiedPhysics::CalculateExternalForcesVersatilCouplingPCISPH(const int i)
{	
	UnifiedParticle &p = particles_[i];

	if (p.type_ == LIQUID_PARTICLE || p.type_ == RIGID_PARTICLE)
	{
		std::vector<int> &indices = neighbor_indices_[i];
		p.force_ = (0.0, 0.0, 0.0);
		const float restVolume = fc_->initialMass  / fc_->fluidRestDensity;

		for (int j = 0; j < indices.size(); j++) {
			UnifiedParticle& neigh = particles_[indices[j]];
			float dist = p.position_.distance(neigh.position_);
			if (dist < fc_->globalSupportRadius) {
				uint dist_lut = dist * (WeightingKernel::lut_size()) / fc_->globalSupportRadius;
				if( dist_lut > WeightingKernel::lut_size())
					dist_lut = WeightingKernel::lut_size();

				// for liquid particles
				if (p.type_ == LIQUID_PARTICLE)
				{
					if (neigh.type_ == LIQUID_PARTICLE)
					{
						// (1) // For each liquid particle, we add viscosity force from other liquid particles	
						// viscosity force & cohesion force
						float kernelVisc = my_kernel_->kernelViscosityLaplacianLut(dist); // symmetric
						// compute artificial viscosity according to MCG03
						// negative symmetry
						vmml::Vector3f v_ij = p.velocity_ - neigh.velocity_;
						p.force_ -= v_ij * restVolume * restVolume * fc_->fluidViscConst * kernelVisc;

						// (2) TODO: surface tension force : cohesion force according to E.q(1) from paper "Versatile Surface Tension and Adhesion for SPH Fluids"
						if (dist >= fc_->globalSupportRadius * 0.25f)	// Avoid Division by Zero Errors
						{
							p.force_ -= CalculateSurfaceTensionForcePCISPHHost(dist_lut, dist, p.position_, neigh.position_);
						}
					} 
					else if (neigh.type_ == FROZEN_PARTICLE || neigh.type_ == RIGID_PARTICLE)
					{
						// For each liquid particle, we add forces from rigid particles & frozen boundary particles
						p.force_ += CalculateBoundaryFluidPressureForceHost(dist_lut, p.position_, neigh.position_, p.predicted_density_, p.previous_correction_pressure_, p.weighted_volume_);

						// TODO: surface tension force : adhesion force according to E.q(6) from paper "Versatile Surface Tension and Adhesion for SPH Fluids"
						if (dist >= fc_->globalSupportRadius * 0.25f)	// Avoid Division by Zero Errors
						{
							p.force_ -= CalculateSurfaceCohesionForcePCISPHHost(dist_lut, dist, p.position_, neigh.position_, neigh.weighted_volume_);
						}
					}
				} 
				else if (p.type_ == RIGID_PARTICLE)
				{
					if (neigh.type_ == RIGID_PARTICLE)	// TODO: add neigh.type == FROZEN_PARTICLE
					{
						// For each rigid particle, we add forces from rigid particles of other rigid bodies
						// "Real-Time Rigid Body Simulation on GPUs" from GPU Gems 3
						const float overlap = 2.0f * fc_->particleRadius - dist;
						// only calculate forces between rigid particles if they belong to different rigid body
						if (overlap > 0.00001f)	// use 0.00001f to smooth numeric error
						{
							// the original one : f = f_spring + f_damping "Real-Time Rigid Body Interaction" (5.15) & (5.16)
							// here r_ij = p_j - p_i  v_ij = v_j - v_i
							p.force_	+= CalculateSpringForceHost(dist, overlap, neigh.position_ - p.position_) + CalculateDampingForceHost(neigh.velocity_ - p.velocity_);
						}
					} 
					else if (neigh.type_ == LIQUID_PARTICLE)
					{
						// For each rigid particle, we add forces from liquid particles using the latter's corrected pressure
						// versatile method E.q(10)
						p.force_ -= CalculateBoundaryFluidPressureForceHost(dist_lut, neigh.position_, p.position_, neigh.predicted_density_, neigh.correction_pressure_, neigh.weighted_volume_);

						if (dist >= fc_->globalSupportRadius * 0.25f)	// Avoid Division by Zero Errors
						{
							p.force_ += CalculateSurfaceCohesionForcePCISPHHost(dist_lut, dist, neigh.position_, p.position_, neigh.weighted_volume_);
						}
					}				
				}			
			}
		}

		// add gravity force to liquid particles
		// Note: for rigid particles, we do not add gravity force, since the gravitational force exerts no torque on rigid body, we delay adding rigid body's 
		// gravitational force to rigid body force & torque calculation
		if (p.type_ == LIQUID_PARTICLE)
		{
			p.force_.y -= fc_->initialMass * fc_->gravityConst;

			// Handle boundary forces exerted on liquid particles
			p.force_ += CalculateBoundaryForcePerLiquidParticleHost(p.position_);

			// Add other external forces in here
		}
		else if (p.type_ == RIGID_PARTICLE)
		{
			// TODO: implement Wall Weight boundary pressure force / versatile boundary particle pressure force & friction force
			// Handle boundary forces exerted on rigid particles
			p.force_ += CalculateBoundaryForcePerRigidParticleHost(p.position_, p.velocity_);

			// Add other external forces in here
		}
		

/*
		// add boundary forces 
		if (fc->addBoundaryForce)
			p.force += BoxBoundaryForce(i);
*/

	}

	p.previous_correction_pressure_ = p.correction_pressure_;
	// init some quantities which are going to be used in the prediction step	
	p.correction_pressure_ = 0.0f;
	p.correction_pressure_force_ = (0.0, 0.0, 0.0);
}

在计算外力的过程中,对于流体粒子,除了计算重力,流体粒子之间的粘性力,还要计算流体粒子之间的内聚力 \mathbf{F}^{cohesion}_{i \leftarrow j} ,公式(16),注意这里并没有实现 \mathbf{F}^{curvature}_{i \leftarrow j}

// (2) TODO: surface tension force : cohesion force according to E.q(1) from paper "Versatile Surface Tension and Adhesion for SPH Fluids"
if (dist >= fc_->globalSupportRadius * 0.25f)	// Avoid Division by Zero Errors
{
	p.force_ -= CalculateSurfaceTensionForcePCISPHHost(dist_lut, dist, p.position_, neigh.position_);
}

vmml::Vector3f UnifiedPhysics::CalculateSurfaceTensionForcePCISPHHost(const uint dist_lut, const float dist, const vmml::Vector3f p_pos, const vmml::Vector3f neigh_pos)
{
	vmml::Vector3f force(0.0f, 0.0f, 0.0f);
	const float splineCoefficient = my_kernel_->lut_spline_surface_tension()[dist_lut];

	force = (p_pos - neigh_pos) * fc_->surface_tension_gamma * fc_->initialMass * fc_->initialMass * splineCoefficient * (1.0f/dist);

	return force;
}

以及流体粒子和固体粒子之间的附着力 \mathbf{F}^{adhesion}_{i \leftarrow k} ,公式(17)

// TODO: surface tension force : adhesion force according to E.q(6) from paper "Versatile Surface Tension and Adhesion for SPH Fluids"
if (dist >= fc_->globalSupportRadius * 0.25f)	// Avoid Division by Zero Errors
{
	p.force_ -= CalculateSurfaceCohesionForcePCISPHHost(dist_lut, dist, p.position_, neigh.position_, neigh.weighted_volume_);
}

vmml::Vector3f UnifiedPhysics::CalculateSurfaceCohesionForcePCISPHHost(const uint dist_lut, const float dist, const vmml::Vector3f p_pos, const vmml::Vector3f neigh_pos, float neigh_weighted_vol)
{
	vmml::Vector3f force(0.0f, 0.0f, 0.0f);
	const float splineCoefficient = my_kernel_->lut_spline_surface_adhesion()[dist_lut];

	force = (p_pos - neigh_pos) * fc_->surface_adhesion_beta * fc_->initialMass * fc_->fluidRestDensity * neigh_weighted_vol * splineCoefficient * (1.0f/dist); 

	return force;
}

FleX实现了上述模型的简化版本,比如表面张力约束用前面介绍的指向流体内部的表面法线 \mathbf{n}_i = \nabla c_s 表示:

C_{tension} = \mathbf{x}_{ij} \cdot \mathbf{n}_i = \cos \theta

Surface Tension Constraint

基本思想跟上面介绍的基于力的模型是一样的,用三个参数进行控制:

	float cohesion;						//!< Control how strongly particles hold each other together, default: 0.025, range [0.0, +inf]
	float surfaceTension;				        //!< Controls how strongly particles attempt to minimize surface area, default: 0.0, range: [0.0, +inf]
        float adhesion;						//!< Controls how strongly particles stick to surfaces they hit, default 0.0, range [0.0, +inf]   

效果如下:

Surface Tensionhttps://www.zhihu.com/video/1047101741941727232adhesionhttps://www.zhihu.com/video/1046870612856938496cohesionhttps://www.zhihu.com/video/1046870666279788544

三、总结:

这篇文章介绍了FleX中流体模拟的理论基础PBF以及流体模拟中特有的表面张力,粘性力,内聚力和附着力等现象的处理方法。相比于PCISPH方法,PBF更稳定,允许大时间步长,更适合应用于游戏物理引擎流体的模拟。另外,基于SPH的方法和基于PBD的方法都非常适合并行化,后续文章讲述GPU并行化的时候再详细介绍。


参考文献:

[1] Macklin, Miles, and Matthias Müller. "Position based fluids." ACM Transactions on Graphics (TOG) 32.4 (2013): 104.

[2] Manual: Concepts

[3] Monaghan, Joseph J. "SPH without a tensile instability." Journal of Computational Physics 159.2 (2000): 290-311.

[4] Fedkiw, Ronald, Jos Stam, and Henrik Wann Jensen. "Visual simulation of smoke." Proceedings of the 28th annual conference on Computer graphics and interactive techniques. ACM, 2001.

[5] Schechter, Hagit, and Robert Bridson. "Ghost SPH for animating water." ACM Transactions on Graphics (TOG) 31.4 (2012): 61.

[6] Akinci, Nadir, et al. "Versatile rigid-fluid coupling for incompressible SPH." ACM Transactions on Graphics (TOG) 31.4 (2012): 62.

[7] Macklin, Miles, et al. "Unified particle physics for real-time applications." ACM Transactions on Graphics (TOG) 33.4 (2014): 153.

[8] Akinci, Nadir, Gizem Akinci, and Matthias Teschner. "Versatile surface tension and adhesion for SPH fluids." ACM Transactions on Graphics (TOG) 32.6 (2013): 182.

[9] Müller, Matthias, David Charypar, and Markus Gross. "Particle-based fluid simulation for interactive applications." Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation. Eurographics Association, 2003.

发布于 2018-11-17

文章被以下专栏收录