Unity物理引擎实战-基于SPH方法的简单水体模拟

franciscochonge · 发表于 2022-2-22 18:39

一个简单的流体效果：

https://www.zhihu.com/video/1471428922994425856
工程：

https://github.com/alen-cell/PhysicsEngine
Navier-Stokes equation方程

基于[1]中的论文的算法实现：
描述一个材料点的受力可以由NS方程给出m理解为f = ma 的ma，一个单位质量乘上加速度。
[2]

来自太极流体力学

关于数学定义：

计算受力

1.pressure

Navier-Stokes equation的第一项便是压强，对于每一个粒子收到的压强不会相同，但是为了使得系统稳定，
将采用一个对称的压强计算

在上一个公式中，计算压强场的梯度之后这个压强的偏差并没有被计算进去。而在计算力之前我们需要计算一个偏移，压强并没有直接被粒子的属性给出，但是可以首先被计算出来：压强的力需要用点所在的位置的密度决定。

ρ0 是静止时水的密度，k是空气含量？根据温度变化而变化

2.Viscosity粘滞力

粘滞力只取决于速度场的偏差而不是绝对的速度，所以计算对称的粘滞力的公式为：

3.表面Tension：

后来再加。。

算法部分：

//ComputeShader kernel定义
private void FindKernels()
{
      clearHashGridKernel = computeShaderSPH.FindKernel(&#34;ClearHashGrid&#34;);
      recalculateHashGridKernel = computeShaderSPH.FindKernel(&#34;RecalculateHashGrid&#34;);
      buildNeighbourListKernel = computeShaderSPH.FindKernel(&#34;BuildNeighbourList&#34;);
      computeDensityPressureKernel = computeShaderSPH.FindKernel(&#34;ComputeDensityPressure&#34;);
      computeForcesKernel = computeShaderSPH.FindKernel(&#34;ComputeForces&#34;);
      integrateKernel = computeShaderSPH.FindKernel(&#34;Integrate&#34;);
      DispatchKernel = computeShaderSPH.FindKernel(&#34;DrawMesh&#34;);

}

核函数

Muller的论文中，给不同的受力不同的核函数，我们使用Wpoly6来做核函数去计算每一个point的物理量(压强和粘滞力除外)：

核函数的梯度和Laplacian算子：

梯度

Laplacian

对于收到的Pressure力，我们用Desbrun的Spiky函数描述，因为这个函数在接近中心的时候的梯度不会突然消失，如果用Wploly6的模型，粒子就会在失去斥力在中心点并集结成Cluster。

在计算压强时

时间积分

所以，我们要计算的便是显示计算受力情况，解出收到的粘滞力和重力之和更新中间速度，在求出压强更新速度，最后用速度更新位置。

Debug

要是看ComputeShader的数值，可以将CommandBuffer的值传递回CPU，
用CommandBuffer.GetData()函数：
//Debug.Log(_argsBuffer);
      //if (material != null && particleMesh != null)
      //{
      Vector3[] data = new Vector3[numberOfParticles];
      _velocitiesBuffer.GetData(data);
      Particle[] dataposition = new Particle[numberOfParticles];
      _particlesBuffer.GetData(dataposition);
      float[] pressureforce = new float[numberOfParticles];
      _pressuresBuffer.GetData(pressureforce);
      float[] densities = new float[numberOfParticles];
      _densitiesBuffer.GetData(densities);


      Vector3[] forces = new Vector3[numberOfParticles];
      _velocitiesBuffer.GetData(forces);

      material.SetFloat(SizeProperty, particleRenderSize);
      material.SetBuffer(ParticlesBufferProperty, _particlesBuffer);

      for (int i = 0; i <Mathf.CeilToInt(numberOfParticles/100); i++)
      {
         //if (pressureforce > 1000|| pressureforce < -1000){
            Debug.Log(i+&#34;number Pressure = &#34; + pressureforce);
         //}
         // Debug.Log(&#34;Velocity = &#34;+ data);
         //Debug.Log(&#34;Position = &#34; + dataposition.position);
         Debug.Log(&#34;Densities = &#34; + densities);
         //Debug.Log(&#34;forces = &#34; + forces);
      }
查看先前定义的数值时可以使用FrameDebugger：

稀疏哈希列表

将将网格的大小设置为粒子的SupportRadius的两倍，只需要遍历周围2X2方向的网格。
也就是每个粒子会拥有2x2x2大小的邻接粒子表。

void GetNearbyKeys(int3 originIndex,float3 position, out int nearbyKeys[8])
{
int3 nearbyBucketIndices[8];
for (int i = 0; i < 8; i++) {
      nearbyBucketIndices = originIndex;
}

if ((originIndex.x + 0.5f) * CellSize <= position.x) {

      nearbyBucketIndices[4].x += 1;
      nearbyBucketIndices[5].x += 1;
      nearbyBucketIndices[6].x += 1;
      nearbyBucketIndices[7].x += 1;

}

else
{
      nearbyBucketIndices[4].x -= 1;
      nearbyBucketIndices[5].x -= 1;
      nearbyBucketIndices[6].x -= 1;
      nearbyBucketIndices[7].x -= 1;
}
if ((originIndex.y + 0.5f) * CellSize <= position.y) {
      nearbyBucketIndices[2].y += 1;
      nearbyBucketIndices[3].y += 1;
      nearbyBucketIndices[6].y += 1;
      nearbyBucketIndices[7].y += 1;
}
else
{
      nearbyBucketIndices[2].y -= 1;
      nearbyBucketIndices[3].y -= 1;
      nearbyBucketIndices[6].y -= 1;
      nearbyBucketIndices[7].y -= 1;
}
if ((originIndex.z + 0.5f) * CellSize <= position.z) {
      nearbyBucketIndices[1].z += 1;
      nearbyBucketIndices[3].z += 1;
      nearbyBucketIndices[5].z += 1;
      nearbyBucketIndices[7].z += 1;
}
else
{
      nearbyBucketIndices[1].z -= 1;
      nearbyBucketIndices[3].z -= 1;
      nearbyBucketIndices[5].z -= 1;
      nearbyBucketIndices[7].z -= 1;
}

for (int j = 0; j < 8; j++) {
      int3 nbcellIndex = nearbyBucketIndices[j];
      if (nbcellIndex.x < 0 || nbcellIndex.x >= Dimensions || nbcellIndex.y < 0 || nbcellIndex.y >= Dimensions || nbcellIndex.z < 0 || nbcellIndex.z >= Dimensions) {
         nearbyKeys[j] = -1;
      }
      else {
         nearbyKeys[j] = Hash(nearbyBucketIndices[j]);
      }
}
}
将每个粒子的临近Particles加入NeighbourList表中：
[numthreads(100, 1, 1)]
void BuildNeighbourList(uint3 id : SV_DispatchThreadID)
{
_neighbourTracker[id.x] = 0;
const int3 cell = GetCell(_particles[id.x].position);
int cells[8];
GetNearbyKeys(cell, _particles[id.x].position, cells);

for (uint j = 0; j < 8; j++)
{
      if (cells[j] == -1) continue; // Grid does not contain cell.
      const uint numberOfParticlesInCell = min(_hashGridTracker[cells[j]], maximumParticlesPerCell); ;
      for (uint index = 0; index < numberOfParticlesInCell; index++)
      {
         const uint potentialNeighbour = _hashGrid[cells[j] * maximumParticlesPerCell + index];
         if (potentialNeighbour == id.x) continue;
         const float3 v = _particles[potentialNeighbour].position - _particles[id.x].position;
         if (dot(v, v) < radius2) // Use squared length (= dot) instead of length for performance.
         {
            _neighbourList[id.x * maximumParticlesPerCell * 8 + _neighbourTracker[id.x]++] = potentialNeighbour;
         }
      }
}
// n. The Neighbouring-list should be n-particles big, each index containing a list of each particles neighbours in radius r.
有了NeighbourList就可以计算压强了，这里的压强我用了Muller论文中的gasConstant来平衡。（大约等于2000f）
因为WCSPH的公式需要小的时间步长。而这里的公式可以达到0.01时间步长不会爆炸。
[numthreads(100,1,1)]//合并
void ComputeDensityPressure(uint3 id :SV_DispatchThreadID)
{
float3 origin = _particles[id.x].position;
float sum = 0.f;
for(int j = 0;j <_neighbourTracker[id.x]; j++)
{
      int neighbourIndex = _neighbourList[id.x*maximumParticlesPerCell*8 +j];
      const float3 diff = origin - _particles[neighbourIndex].position;
      float distanceSquared = dot(diff,diff);
      sum += Wpoly6Kernel(distanceSquared);
}
//computeDensity

//TodO presure
_densities[id.x] = sum*mass + 0.000001f;
_densities[id.x] = max(restDensity, _densities[id.x]);
_pressures[id.x] =  gasConstant * (_densities[id.x] - restDensity);

}Tips：（若是不维护一个固定的表格，也可以用一个RadixSort实现，但是还不知道排序的时间和这种方式哪个更省，先都探索一下）参考：https://developer.download.nvidia.cn/assets/cuda/files/particles.pdf

进一步解决的问题：

1.粒子在边界收到的压强会小于应该有的压强，导致粒子粘滞在边界。这个需要通过虚粒子的方式解决掉。
2.碰撞。碰撞有两种方式：SDF,固体粒子。而这两种方法在后面都会加入系统之中。因为两种方法其实密不可分，一种重建粒子表面的方法便是建立在SDF之上的。
3.流体的不可压缩性。将使用PBD计算来解决。PositioBasedDynamics。
4.流体和固体的双向耦合。统一用粒子方式解决。
性能：

20000-40000个粒子稳定在160帧左右：
后期加上屏幕渲染方式，固体粒子双向耦合希望可以在60-100 帧上下

资料：

下一篇：基于PBD的流体模拟和刚体耦合
参考

^1https://matthias-research.github.io/pages/publications/sca03.pdf

^2https://www.bilibili.com/video/BV1mi4y1o7wz?from=search&seid=16843922597556751206&spm_id_from=333.337.0.0

闲鱼技术01 · 发表于 2022-2-22 18:40

我什么才能像你这般优秀[捂脸]

Ylisar · 发表于 2022-2-22 18:41

我也有很多不懂得，刚开始学习[调皮]

		自动登录	找回密码
密码			立即注册

Unity物理引擎实战-基于SPH方法的简单水体模拟

本帖子中包含更多资源