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

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

一个简单的流体效果：




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

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

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




来自太极 流体力学



关于数学定义：




计算受力

1.pressure



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


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



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

2.Viscosity粘滞力

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


3.表面Tension：

后来再加。。


算法部分：


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

    }

核函数

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;
      _velocitiesBuffer.GetData(data);
      Particle[] dataposition = new Particle;
      _particlesBuffer.GetData(dataposition);
      float[] pressureforce = new float;
      _pressuresBuffer.GetData(pressureforce);
      float[] densities = new float;
      _densitiesBuffer.GetData(densities);
      
      
      Vector3[] forces = new Vector3;
      _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+"number Pressure = " + pressureforce);
            //}
            // Debug.Log("Velocity = "+ data);
            //Debug.Log("Position = " + dataposition.position);
            Debug.Log("Densities = " + densities);
            //Debug.Log("forces = " + forces);
      }
查看先前定义的数值时可以使用FrameDebugger：



稀疏哈希列表

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




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


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

    }

    else
    {
      nearbyBucketIndices.x -= 1;
      nearbyBucketIndices.x -= 1;
      nearbyBucketIndices.x -= 1;
      nearbyBucketIndices.x -= 1;
    }
    if ((originIndex.y + 0.5f) * CellSize <= position.y) {
      nearbyBucketIndices.y += 1;
      nearbyBucketIndices.y += 1;
      nearbyBucketIndices.y += 1;
      nearbyBucketIndices.y += 1;
    }
    else
    {
      nearbyBucketIndices.y -= 1;
      nearbyBucketIndices.y -= 1;
      nearbyBucketIndices.y -= 1;
      nearbyBucketIndices.y -= 1;
    }
    if ((originIndex.z + 0.5f) * CellSize <= position.z) {
      nearbyBucketIndices.z += 1;
      nearbyBucketIndices.z += 1;
      nearbyBucketIndices.z += 1;
      nearbyBucketIndices.z += 1;
    }
    else
    {
      nearbyBucketIndices.z -= 1;
      nearbyBucketIndices.z -= 1;
      nearbyBucketIndices.z -= 1;
      nearbyBucketIndices.z -= 1;
    }

    for (int j = 0; j < 8; j++) {
      int3 nbcellIndex = nearbyBucketIndices;
      if (nbcellIndex.x < 0 || nbcellIndex.x >= Dimensions || nbcellIndex.y < 0 || nbcellIndex.y >= Dimensions || nbcellIndex.z < 0 || nbcellIndex.z >= Dimensions) {
            nearbyKeys = -1;
      }
      else {
            nearbyKeys = Hash(nearbyBucketIndices);
      }
    }
}
将每个粒子的临近Particles加入NeighbourList表中：

void BuildNeighbourList(uint3 id : SV_DispatchThreadID)
{
    _neighbourTracker = 0;
    const int3 cell = GetCell(_particles.position);
    int cells;
    GetNearbyKeys(cell, _particles.position, cells);

    for (uint j = 0; j < 8; j++)
    {
      if (cells == -1) continue; // Grid does not contain cell.
      const uint numberOfParticlesInCell = min(_hashGridTracker], maximumParticlesPerCell); ;
      for (uint index = 0; index < numberOfParticlesInCell; index++)
      {
            const uint potentialNeighbour = _hashGrid * maximumParticlesPerCell + index];
            if (potentialNeighbour == id.x) continue;
            const float3 v = _particles.position - _particles.position;
            if (dot(v, v) < radius2) // Use squared length (= dot) instead of length for performance.
            {
                _neighbourList++] = 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时间步长不会爆炸。
//合并
void ComputeDensityPressure(uint3 id :SV_DispatchThreadID)
{   
    float3 origin = _particles.position;
    float sum = 0.f;
    for(int j = 0;j <_neighbourTracker; j++)
{      
      int neighbourIndex = _neighbourList;
      const float3 diff = origin - _particles.position;
      float distanceSquared = dot(diff,diff);
      sum += Wpoly6Kernel(distanceSquared);
}
//computeDensity

    //TodO presure
_densities = sum*mass + 0.000001f;
_densities = max(restDensity, _densities);
_pressures =gasConstant * (_densities - 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

我也有很多不懂得，刚开始学习[调皮]
页: [1]
查看完整版本: Unity物理引擎实战-基于SPH方法的简单水体模拟