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

franciscochonge

一个简单的流体效果：




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("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[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+"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[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

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

Ylisar

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