技术美术｜游戏中的流体模拟（Fluid Simulation）

技术美术｜游戏中的流体模拟（Fluid Simulation）

  • 作者：admin
  • /
  • 时间：2023年12月06日
  • /
  • 浏览：834 次
  • /
  • 分类：厚积薄发

【USparkle专栏】如果你深怀绝技，爱“搞点研究”，乐于分享也博采众长，我们期待你的加入，让智慧的火花碰撞交织，让知识的传递生生不息！

一、闲聊

最近一直在研究流体模拟，很神奇的一个东西，好在网上有很多参考资料，研究过程不算太困难。分享下最近一段时间的学习心得。

二、效果演示

三、算法原理

游戏领域实现流体模拟的几种常见方式有：

  • 基于网格的方法：在网格上模拟，每个格子都有自己的数据（速度、密度、颜色、温度等），逐帧更新格子内数据。这种方法的优点是方便多线程实现，渲染也很方便。缺点是计算过程中需要对参数做估算，容易产生误差。
  • 基于粒子的方法：将流体具象化为很多个粒子，每个粒子都有自己数据（速度、颜色、温度），逐帧更新粒子的位置。这种方法的优点是误差小，能表现出更多的流体细节。缺点是不利于多线程实现，渲染也比较麻烦。

这篇文章采用的是基于网格的方法，流体有很多种类（气体、水、岩浆、蜂蜜等），不同流体使用的算法各有差异，这篇文章讨论的是气体流体模拟。

在流体模拟中，有两个主要计算过程，压缩解算和流动。

压缩解算（Project）
压缩性是流体的基本属性之一，正常环境下，大多数流体都很难被压缩，向流体施加很大的力，而流体的体积变化却很小。

压缩解算的目的，就是要模拟流体很难被压缩的特点，假设我们在一个8x8的网格上做流体模拟：

先在格子边框上创建辅助点（Staggered Point），水平方向辅助点为黄色，垂直方向辅助点为绿色：

拿中间几个格子举例，每个格子都有自己的速度：

将格子的速度拆分到周围4个辅助点上，水平速度存入黄色点，垂直速度存入绿色点：


单个格子拆分前


单个格子拆分后


整体拆分前


整体拆分后

然后根据格子周围4个辅助点的速度，对格子做压缩解算：

上图的格子有三个方向在流入，一个方向在流出，流入量大于流出量，要使流体不被压缩，流入量和流出量必须相等。

先计算净流入、流出量（Divergence）：

float divergence = leftPointSpeed - rightPointSpeed + downPointSpeed - upPointSpeed;

将其均分后修改辅助点速度：

divergence /= 4;
leftPointSpeed += -divergence;
rightPointSpeed += divergence;
downPointSpeed += -divergence;
upPointSpeed += divergence;


修改辅助点速度后

这样就保证了这个格子流入和流出量相等。

再看一个流体遇到障碍物的例子：

格子的右边是一面墙，所以右边黄色点的速度始终为0，压缩解算的公式变为：

float divergence = leftPointSpeed + downPointSpeed - upPointSpeed;

divergence /= 3;
leftPointSpeed += -divergence;
downPointSpeed += -divergence;
upPointSpeed += divergence;


修改辅助点速度后

我们可以为每个辅助点附加一个Scaler，障碍物的Scaler为0，非障碍物的Scaler为1，这样一来，有无障碍物都可以使用同一个公式：

int counter = leftPointScaler + rightPointScaler + downPointScaler + upPointScaler; 

float divergence = leftPointSpeed * leftPointScaler 
                   - rightPointSpeed * rightPointScaler 
                   + downPointSpeed * downPointScaler 
                   - upPointSpeed * upPointScaler;

divergence /= counter;

leftPointSpeed += -divergence * leftPointScaler;
rightPointSpeed += divergence * rightPointScaler;
downPointSpeed += -divergence * downPointScaler;
upPointSpeed += divergence * upPointScaler;

最后，计算辅助点速度的平均值，更新格子速度：

float uSpeed = (leftPointSpeed + rightPointSpeed) / 2;
float vSpeed = (downPointSpeed + upPointSpeed) / 2;

cellData.speed = float2(uSpeed, vSpeed);

一个格子速度发生变化，其临近格子的流入、流出量也会改变，这里我们需要迭代多次去逼近正确解：

for(int i = 0; i < iteration; i++) {
    for(int j = 0; j < cellNumber; j++) {
        //对格子做计算，修改辅助点速度
    }

    //更新格子速度
}

这样，压缩解算就完成了！

流动（Advect）
更新完格子的速度后，就可以移动格子内的数据了。最直观的做法是，根据格子的速度，计算出它移动到了哪个位置，然后把它的数据（密度，速度等）加入到新格子中。

这种做法最直观，很好理解，但存在一个问题，可能会有多个格子移动到了同一个位置：

在多线程计算时，对新格子数据读写会出现冲突。要解决这个问题，通常采用的方法是逆向过来，先估算格子的速度，反过来找到它在移动前的位置，用移动前位置周围几个格子内的数据做插值，更新自己。

估算速度可以用格子和其周围8个格子速度的平均值：

也可以用周围12个辅助点速度的加权平均值：

扩散（Diffuse）
流体还有扩散特性，高浓度区域会主动扩散到低浓度区域，直至所有格子的浓度相等。比方说向水杯里滴入一滴墨水，墨水会逐渐扩散开，直至整杯水均匀变黑。不过这一步并不是必需的，我在实际尝试中发现，加入扩散后效果反而没那么好看了。

四、Unity内具体实现过程

我使用的Unity版本是2021.3，URP管线。流体模拟的计算量比较大，我使用的ComputeShader做计算。

主要流程

private void OnEnable() {
    //定义并初始化数据结构
}

private void Update() {
    //向指定格子输入流体

    //将格子速度拆分到辅助点上

    //压缩解算，修改辅助点速度

    //更新格子速度

    //格子数据流动

    //衰减
}

数据结构
CellData为单个格子的数据结构，UStaggeredPoint代表水平方向的辅助点，VStaggeredPoint代表垂直方向的辅助点。

int2 _Resolution;

struct CellData {
    int2 coord;
    float density;
    float2 velocity;
    float4 color;
    int2 leftStaggeredPointCoord;
    int2 rightStaggeredPointCoord;
    int2 upStaggeredPointCoord;
    int2 downStaggeredPointCoord;
    int leftStaggeredPointIndex;
    int rightStaggeredPointIndex;
    int upStaggeredPointIndex;
    int downStaggeredPointIndex;
    int leftStaggeredPointSummaryIndex;
    int rightStaggeredPointSummaryIndex;
    int upStaggeredPointSummaryIndex;
    int downStaggeredPointSummaryIndex;
};

struct StaggeredPoint {
    int2 coord;
    float scaler;
    float velocity;
    int summaryNumber;
};

int CellCoordToIndex(int2 coord) {
    return coord.x + coord.y * _Resolution.x;
}

int UStaggeredPointCoordToIndex(int2 coord) {
    return coord.x + coord.y * (_Resolution.x + 1);
}

int VStaggeredPointCoordToIndex(int2 coord) {
    return coord.x + coord.y * _Resolution.x;
}

#endif

初始化
初始化格子数据：

#pragma kernel CSMain

#include "Assets/FluidSimulationLibrary.hlsl"

RWStructuredBuffer<CellData> _CellDatas;
RWStructuredBuffer<StaggeredPoint> _UStaggeredPoints;
RWStructuredBuffer<StaggeredPoint> _VStaggeredPoints;

[numthreads(1,1,1)]
void CSMain (uint3 id : SV_DispatchThreadID)
{
uint index = id.x + id.y * _Resolution.x;
if(index >= _CellDatas.Length) {
return;
}

CellData cellData = _CellDatas[index];

cellData.coord = id.xy;
cellData.density = 0;
cellData.velocity = 0;

        int2 leftStaggeredPointCoord = cellData.coord + int2(0, 0);
        int2 rightStaggeredPointCoord = cellData.coord + int2(1, 0);
        int2 upStaggeredPointCoord = cellData.coord + int2(0, 1);
        int2 downStaggeredPointCoord = cellData.coord + int2(0, 0);

cellData.leftStaggeredPointCoord = leftStaggeredPointCoord;
cellData.rightStaggeredPointCoord = rightStaggeredPointCoord;
cellData.upStaggeredPointCoord = upStaggeredPointCoord;
cellData.downStaggeredPointCoord = downStaggeredPointCoord;

        int leftStaggeredPointIndex = UStaggeredPointCoordToIndex(leftStaggeredPointCoord);
        int rightStaggeredPointIndex = UStaggeredPointCoordToIndex(rightStaggeredPointCoord);
        int upStaggeredPointIndex = VStaggeredPointCoordToIndex(upStaggeredPointCoord);
        int downStaggeredPointIndex = VStaggeredPointCoordToIndex(downStaggeredPointCoord);

cellData.leftStaggeredPointIndex = leftStaggeredPointIndex;
cellData.rightStaggeredPointIndex = rightStaggeredPointIndex;
cellData.upStaggeredPointIndex = upStaggeredPointIndex;
cellData.downStaggeredPointIndex = downStaggeredPointIndex;

        int leftStaggeredPointSummaryNumber;
        int rightStaggeredPointSummaryNumber;
        int upStaggeredPointSummaryNumber;
        int downStaggeredPointSummaryNumber;

InterlockedAdd(_UStaggeredPoints[leftStaggeredPointIndex].summaryNumber, 1, leftStaggeredPointSummaryNumber);
InterlockedAdd(_UStaggeredPoints[rightStaggeredPointIndex].summaryNumber, 1, rightStaggeredPointSummaryNumber);
InterlockedAdd(_VStaggeredPoints[upStaggeredPointIndex].summaryNumber, 1, upStaggeredPointSummaryNumber);
InterlockedAdd(_VStaggeredPoints[downStaggeredPointIndex].summaryNumber, 1, downStaggeredPointSummaryNumber);

cellData.leftStaggeredPointSummaryIndex = leftStaggeredPointIndex * 2 + leftStaggeredPointSummaryNumber;
cellData.rightStaggeredPointSummaryIndex = rightStaggeredPointIndex * 2 + rightStaggeredPointSummaryNumber;
cellData.upStaggeredPointSummaryIndex = (upStaggeredPointIndex + _UStaggeredPoints.Length) * 2 + upStaggeredPointSummaryNumber;
cellData.downStaggeredPointSummaryIndex = (downStaggeredPointIndex + _UStaggeredPoints.Length) * 2 + downStaggeredPointSummaryNumber;

_CellDatas[index] = cellData;
}

初始化辅助点数据，C#会调用ComputeShader两次，分别初始化水平和垂直方向的辅助点：

#pragma kernel CSMain

#include "Assets/FluidSimulationLibrary.hlsl"

RWStructuredBuffer<StaggeredPoint> _StaggeredPoints;

int _ColumnNumber;
int _WallThickness;

[numthreads(1,1,1)]
void CSMain (uint3 id : SV_DispatchThreadID)
{
    uint index = id.x + id.y * _ColumnNumber;
    if(index >= _StaggeredPoints.Length) {
        return;
    }

    StaggeredPoint staggeredPoint = _StaggeredPoints[index];
    staggeredPoint.coord = id.xy;
    staggeredPoint.scaler = 1;
    staggeredPoint.velocity = 0;
    staggeredPoint.summaryNumber = 0;

    if(_ColumnNumber == _Resolution.x) {
        if(staggeredPoint.coord.y < _WallThickness) {
    staggeredPoint.scaler = 0;
}
else if(staggeredPoint.coord.y > _Resolution.y - _WallThickness) {
    staggeredPoint.scaler = 0;
}
    }
    else {
if(staggeredPoint.coord.x < _WallThickness) {
    staggeredPoint.scaler = 0;
}
else if(staggeredPoint.coord.x > _Resolution.x - _WallThickness) {
    staggeredPoint.scaler = 0;
}
    } 

    _StaggeredPoints[index] = staggeredPoint;
}

输入
当按下鼠标左键时，通过C#将输入信息传入ComputeShader，找到鼠标周围的格子，修改数据：

struct InjectData {
    int2 center;
    float radius;
    float density;
    float2 velocity;
    float4 color;
};

#pragma kernel CSMain

#include "Assets/FluidSimulationLibrary.hlsl"

RWStructuredBuffer<CellData> _CellDatas;
StructuredBuffer<InjectData> _InjectDatas;
StructuredBuffer<StaggeredPoint> _UStaggeredPoints;
StructuredBuffer<StaggeredPoint> _VStaggeredPoints;

[numthreads(1,1,1)]
void CSMain (uint3 id : SV_DispatchThreadID)
{
    uint index = id.x + id.y * _Resolution.x;
    if(index >= _CellDatas.Length) {
        return;
    }

    CellData cellData = _CellDatas[index];
    for(uint iii = 0; iii < _InjectDatas.Length; iii++) {
        InjectData injectData = _InjectDatas[iii];

float dist = distance(cellData.coord, injectData.center);
float t = 1 - saturate(dist / injectData.radius);

if(t > 0) {
    StaggeredPoint leftStaggeredPoint = _UStaggeredPoints[cellData.leftStaggeredPointIndex];
    StaggeredPoint rightStaggeredPoint = _UStaggeredPoints[cellData.rightStaggeredPointIndex];
    StaggeredPoint upStaggeredPoint = _VStaggeredPoints[cellData.upStaggeredPointIndex];
    StaggeredPoint downStaggeredPoint = _VStaggeredPoints[cellData.downStaggeredPointIndex];

    if(leftStaggeredPoint.scaler == 0 || rightStaggeredPoint.scaler == 0 || upStaggeredPoint.scaler == 0 || downStaggeredPoint.scaler == 0) {
        continue;
    }

    cellData.density += injectData.density * t;
    cellData.velocity += injectData.velocity * t;
    cellData.color += injectData.color * t;
}
    }

    _CellDatas[index] = cellData;
}

速度拆分
将格子的速度拆分到辅助点上，先累加，再求平均：

#pragma kernel CSMain

#include "Assets/FluidSimulationLibrary.hlsl"

StructuredBuffer<CellData> _CellDatas;
RWStructuredBuffer<float> _SummaryDatas;

[numthreads(256,1,1)]
void CSMain (uint3 id : SV_DispatchThreadID)
{
    uint index = id.x;
    if(index >= _CellDatas.Length) {
        return;
    }

    CellData cellData = _CellDatas[index];

    float2 velocity = cellData.velocity;

    if(velocity.x > 0) {
        _SummaryDatas[cellData.rightStaggeredPointSummaryIndex] = velocity.x;
    }
    else if(velocity.x < 0) {
        _SummaryDatas[cellData.leftStaggeredPointSummaryIndex] = velocity.x;
    }

    if(velocity.y > 0) {
        _SummaryDatas[cellData.upStaggeredPointSummaryIndex] = velocity.y;
    }
    else if(velocity.y < 0) {
        _SummaryDatas[cellData.downStaggeredPointSummaryIndex] = velocity.y;
    }
}

#pragma kernel CSMain

#include "Assets/FluidSimulationLibrary.hlsl"

RWStructuredBuffer<StaggeredPoint> _UStaggeredPoints;
RWStructuredBuffer<StaggeredPoint> _VStaggeredPoints;
RWStructuredBuffer<float> _SummaryDatas;

float AverageVelocity(StaggeredPoint staggeredPoint, int index) {
    int counter = 0;
    float velocity = 0;

    int summaryIndex0 = index * 2;
    int summaryIndex1 = index * 2 + 1;

    if(_SummaryDatas[summaryIndex0] != 0) {
        counter += 1;
        velocity += _SummaryDatas[summaryIndex0];
        _SummaryDatas[summaryIndex0] = 0;
    }

    if(_SummaryDatas[summaryIndex1] != 0) {
        counter += 1;
        velocity += _SummaryDatas[summaryIndex1];
       _SummaryDatas[summaryIndex1] = 0;
    }

    if(counter == 0) {
        return 0;
    }
    else {
        return velocity / counter;
    }
}

[numthreads(256,1,1)]
void CSMain (uint3 id : SV_DispatchThreadID)
{
    uint index = id.x;
    if(index >= _UStaggeredPoints.Length + _VStaggeredPoints.Length) {
return;
}

    StaggeredPoint staggeredPoint;

    if(index >= _UStaggeredPoints.Length) {
        staggeredPoint = _VStaggeredPoints[index % _UStaggeredPoints.Length];
        staggeredPoint.velocity = AverageVelocity(staggeredPoint, index) * staggeredPoint.scaler;
        _VStaggeredPoints[index % _UStaggeredPoints.Length] = staggeredPoint;
    }
    else {
        staggeredPoint = _UStaggeredPoints[index];
        staggeredPoint.velocity = AverageVelocity(staggeredPoint, index) * staggeredPoint.scaler;
        _UStaggeredPoints[index] = staggeredPoint;
    }
}

压缩解算
根据净流入、流出量修改辅助点速度，先累加，再求平均：

#pragma kernel CSMain

#include "Assets/FluidSimulationLibrary.hlsl"

StructuredBuffer<CellData> _CellDatas;
StructuredBuffer<StaggeredPoint> _UStaggeredPoints;
StructuredBuffer<StaggeredPoint> _VStaggeredPoints;
RWStructuredBuffer<float> _SummaryDatas;

[numthreads(256,1,1)]
void CSMain (uint3 id : SV_DispatchThreadID)
{
    uint index = id.x;
    if(index >= _CellDatas.Length) {
        return;
    }

    CellData cellData = _CellDatas[index];

    StaggeredPoint leftStaggeredPoint = _UStaggeredPoints[cellData.leftStaggeredPointIndex];
    StaggeredPoint rightStaggeredPoint = _UStaggeredPoints[cellData.rightStaggeredPointIndex];
    StaggeredPoint upStaggeredPoint = _VStaggeredPoints[cellData.upStaggeredPointIndex];
    StaggeredPoint downStaggeredPoint = _VStaggeredPoints[cellData.downStaggeredPointIndex];

    int leftScaler = leftStaggeredPoint.scaler;
    int rightScaler = rightStaggeredPoint.scaler;
    int upScaler = upStaggeredPoint.scaler;
    int downScaler = downStaggeredPoint.scaler;

    int counter = (leftScaler + rightScaler + upScaler + downScaler);

    if(counter == 0) {
        return;
    }

    float divergence = (leftStaggeredPoint.velocity 
        - rightStaggeredPoint.velocity 
        - upStaggeredPoint.velocity 
        + downStaggeredPoint.velocity) 
        / counter;

    _SummaryDatas[cellData.leftStaggeredPointSummaryIndex] = -divergence * leftScaler;
    _SummaryDatas[cellData.rightStaggeredPointSummaryIndex] = divergence * rightScaler;
    _SummaryDatas[cellData.upStaggeredPointSummaryIndex] = divergence * upScaler;
    _SummaryDatas[cellData.downStaggeredPointSummaryIndex] = -divergence * downScaler;
}

#pragma kernel CSMain

#include "Assets/FluidSimulationLibrary.hlsl"

RWStructuredBuffer<StaggeredPoint> _UStaggeredPoints;
RWStructuredBuffer<StaggeredPoint> _VStaggeredPoints;
RWStructuredBuffer<float> _SummaryDatas;

float AverageVelocity(StaggeredPoint staggeredPoint, int index) {
    int counter = 0;
    float velocity = 0;

    int summaryIndex0 = index * 2;
    int summaryIndex1 = index * 2 + 1;

    if(_SummaryDatas[summaryIndex0] != 0) {
        counter += 1;
        velocity += _SummaryDatas[summaryIndex0];
        _SummaryDatas[summaryIndex0] = 0;
    }

    if(_SummaryDatas[summaryIndex1] != 0) {
        counter += 1;
        velocity += _SummaryDatas[summaryIndex1];
       _SummaryDatas[summaryIndex1] = 0;
    }

    if(counter == 0) {
        return staggeredPoint.velocity;
    }
    else {
        return staggeredPoint.velocity + velocity / counter;
    }
}

[numthreads(256,1,1)]
void CSMain (uint3 id : SV_DispatchThreadID)
{
    uint index = id.x;
    if(index >= _UStaggeredPoints.Length + _VStaggeredPoints.Length) {
return;
}

    StaggeredPoint staggeredPoint;

    if(index >= _UStaggeredPoints.Length) {
        staggeredPoint = _VStaggeredPoints[index % _UStaggeredPoints.Length];
        staggeredPoint.velocity = AverageVelocity(staggeredPoint, index);
        _VStaggeredPoints[index % _UStaggeredPoints.Length] = staggeredPoint;
    }
    else {
        staggeredPoint = _UStaggeredPoints[index];
        staggeredPoint.velocity = AverageVelocity(staggeredPoint, index);
        _UStaggeredPoints[index] = staggeredPoint;
    }
}

更新格子速度

#pragma kernel CSMain

#include "Assets/FluidSimulationLibrary.hlsl"

RWStructuredBuffer<CellData> _CellDatas;
StructuredBuffer<StaggeredPoint> _UStaggeredPoints;
StructuredBuffer<StaggeredPoint> _VStaggeredPoints;

[numthreads(256,1,1)]
void CSMain (uint3 id : SV_DispatchThreadID)
{
    uint index = id.x;
    if(index >= _CellDatas.Length) {
        return;
    }

    CellData cellData = _CellDatas[index];

    StaggeredPoint leftStaggeredPoint = _UStaggeredPoints[cellData.leftStaggeredPointIndex];
    StaggeredPoint rightStaggeredPoint = _UStaggeredPoints[cellData.rightStaggeredPointIndex];
    StaggeredPoint upStaggeredPoint = _VStaggeredPoints[cellData.upStaggeredPointIndex];
    StaggeredPoint downStaggeredPoint = _VStaggeredPoints[cellData.downStaggeredPointIndex];

    cellData.velocity.x = (leftStaggeredPoint.velocity + rightStaggeredPoint.velocity) / 2;
    cellData.velocity.y = (upStaggeredPoint.velocity + downStaggeredPoint.velocity) / 2;

    _CellDatas[index] = cellData;
}

流动
先估算格子速度，用格子周围12个辅助点的加权平均值。算出格子在流动前的位置，对流动前位置临近4个格子内的数据做插值，更新自己：

struct AdvectData {
    float density;
    float2 velocity;
    float4 color;
};

#pragma kernel CSMain

#include "Assets/FluidSimulationLibrary.hlsl"

StructuredBuffer<CellData> _CellDatas;
RWStructuredBuffer<AdvectData> _AdvectDatas;
StructuredBuffer<StaggeredPoint> _UStaggeredPoints;
StructuredBuffer<StaggeredPoint> _VStaggeredPoints;

[numthreads(256,1,1)]
void CSMain (uint3 id : SV_DispatchThreadID)
{
    uint index = id.x;
    if(index >= _CellDatas.Length) {
        return;
    }

    CellData cellData = _CellDatas[index];

    float uVelocity = 0;
    float vVelocity = 0;
    int uCounter = 0;
    int vCounter = 0;

    int2 pointCoord = cellData.leftStaggeredPointCoord + int2(0, 0);
    if(pointCoord.x >= 0 && pointCoord.x <= _Resolution.x && pointCoord.y >= 0 && pointCoord.y < _Resolution.y) {
        int pointIndex = UStaggeredPointCoordToIndex(pointCoord);
StaggeredPoint staggeredPoint = _UStaggeredPoints[pointIndex];
uVelocity += staggeredPoint.velocity * staggeredPoint.scaler * 2;
uCounter += staggeredPoint.scaler * 2;
    }

    pointCoord = cellData.leftStaggeredPointCoord + int2(0, 1);
    if(pointCoord.x >= 0 && pointCoord.x <= _Resolution.x && pointCoord.y >= 0 && pointCoord.y < _Resolution.y) {
        int pointIndex = UStaggeredPointCoordToIndex(pointCoord);
StaggeredPoint staggeredPoint = _UStaggeredPoints[pointIndex];
uVelocity += staggeredPoint.velocity * staggeredPoint.scaler;
uCounter += staggeredPoint.scaler;
    }

    pointCoord = cellData.leftStaggeredPointCoord + int2(0, -1);
    if(pointCoord.x >= 0 && pointCoord.x <= _Resolution.x && pointCoord.y >= 0 && pointCoord.y < _Resolution.y) {
        int pointIndex = UStaggeredPointCoordToIndex(pointCoord);
StaggeredPoint staggeredPoint = _UStaggeredPoints[pointIndex];
uVelocity += staggeredPoint.velocity * staggeredPoint.scaler;
uCounter += staggeredPoint.scaler;
    }

    pointCoord = cellData.rightStaggeredPointCoord + int2(0, 0);
    if(pointCoord.x >= 0 && pointCoord.x <= _Resolution.x && pointCoord.y >= 0 && pointCoord.y < _Resolution.y) {
        int pointIndex = UStaggeredPointCoordToIndex(pointCoord);
StaggeredPoint staggeredPoint = _UStaggeredPoints[pointIndex];
uVelocity += staggeredPoint.velocity * staggeredPoint.scaler * 2;
uCounter += staggeredPoint.scaler * 2;
    }

    pointCoord = cellData.rightStaggeredPointCoord + int2(0, 1);
    if(pointCoord.x >= 0 && pointCoord.x <= _Resolution.x && pointCoord.y >= 0 && pointCoord.y < _Resolution.y) {
int pointIndex = UStaggeredPointCoordToIndex(pointCoord);
StaggeredPoint staggeredPoint = _UStaggeredPoints[pointIndex];
uVelocity += staggeredPoint.velocity * staggeredPoint.scaler;
uCounter += staggeredPoint.scaler;
    }

    pointCoord = cellData.rightStaggeredPointCoord + int2(0, -1);
    if(pointCoord.x >= 0 && pointCoord.x <= _Resolution.x && pointCoord.y >= 0 && pointCoord.y < _Resolution.y) {
int pointIndex = UStaggeredPointCoordToIndex(pointCoord);
StaggeredPoint staggeredPoint = _UStaggeredPoints[pointIndex];
uVelocity += staggeredPoint.velocity * staggeredPoint.scaler;
uCounter += staggeredPoint.scaler;
    }

    pointCoord = cellData.upStaggeredPointCoord + int2(0, 0);
    if(pointCoord.x >= 0 && pointCoord.x < _Resolution.x && pointCoord.y >= 0 && pointCoord.y <= _Resolution.y) {
int pointIndex = VStaggeredPointCoordToIndex(pointCoord);
StaggeredPoint staggeredPoint = _VStaggeredPoints[pointIndex];
vVelocity += staggeredPoint.velocity * staggeredPoint.scaler * 2;
vCounter += staggeredPoint.scaler * 2;
    }

    pointCoord = cellData.upStaggeredPointCoord + int2(1, 0);
    if(pointCoord.x >= 0 && pointCoord.x < _Resolution.x && pointCoord.y >= 0 && pointCoord.y <= _Resolution.y) {
int pointIndex = VStaggeredPointCoordToIndex(pointCoord);
StaggeredPoint staggeredPoint = _VStaggeredPoints[pointIndex];
vVelocity += staggeredPoint.velocity * staggeredPoint.scaler;
vCounter += staggeredPoint.scaler;
    }

    pointCoord = cellData.upStaggeredPointCoord + int2(-1, 0);
    if(pointCoord.x >= 0 && pointCoord.x < _Resolution.x && pointCoord.y >= 0 && pointCoord.y <= _Resolution.y) {
int pointIndex = VStaggeredPointCoordToIndex(pointCoord);
StaggeredPoint staggeredPoint = _VStaggeredPoints[pointIndex];
vVelocity += staggeredPoint.velocity * staggeredPoint.scaler;
vCounter += staggeredPoint.scaler;
    }

    pointCoord = cellData.downStaggeredPointCoord + int2(0, 0);
    if(pointCoord.x >= 0 && pointCoord.x < _Resolution.x && pointCoord.y >= 0 && pointCoord.y <= _Resolution.y) {
int pointIndex = VStaggeredPointCoordToIndex(pointCoord);
StaggeredPoint staggeredPoint = _VStaggeredPoints[pointIndex];
vVelocity += staggeredPoint.velocity * staggeredPoint.scaler * 2;
vCounter += staggeredPoint.scaler * 2;
    }

    pointCoord = cellData.downStaggeredPointCoord + int2(1, 0);
    if(pointCoord.x >= 0 && pointCoord.x < _Resolution.x && pointCoord.y >= 0 && pointCoord.y <= _Resolution.y) {
int pointIndex = VStaggeredPointCoordToIndex(pointCoord);
StaggeredPoint staggeredPoint = _VStaggeredPoints[pointIndex];
vVelocity += staggeredPoint.velocity * staggeredPoint.scaler;
vCounter += staggeredPoint.scaler;
    }

    pointCoord = cellData.downStaggeredPointCoord + int2(-1, 0);
    if(pointCoord.x >= 0 && pointCoord.x < _Resolution.x && pointCoord.y >= 0 && pointCoord.y <= _Resolution.y) {
int pointIndex = VStaggeredPointCoordToIndex(pointCoord);
StaggeredPoint staggeredPoint = _VStaggeredPoints[pointIndex];
vVelocity += staggeredPoint.velocity * staggeredPoint.scaler;
vCounter += staggeredPoint.scaler;
    }

    if(uCounter == 0) {
        uVelocity = 0;
    }
    else {
        uVelocity /= uCounter;
    }

    if(vCounter == 0) {
        vVelocity = 0;
    }
    else {
vVelocity /= vCounter;
    }



    float ut;
    float vt;
    int leftX;
    int rightX;
    int upY;
    int downY;

    float udist = -uVelocity;
    if(udist > 0) {
ut = frac(udist);
        leftX = cellData.coord.x + floor(udist);
        rightX = leftX + 1;
leftX = min(leftX, _Resolution.x - 1);
        rightX = min(rightX, _Resolution.x - 1);
    }
    else {
        udist = abs(udist);
        ut = 1 - frac(udist);
        leftX = cellData.coord.x - floor(udist) - 1;
        rightX = leftX + 1;
        leftX = max(leftX, 0);
        rightX = max(rightX, 0);
    }

    float vdist = -vVelocity;
    if(vdist > 0) {
vt = frac(vdist);
        downY = cellData.coord.y + floor(vdist);
        upY = downY + 1;
downY = min(downY, _Resolution.y - 1);
        upY = min(upY, _Resolution.y - 1);
    }
    else {
        vdist = abs(vdist);
        vt = 1 - frac(vdist);
        downY = cellData.coord.y - floor(vdist) - 1;
        upY = downY + 1;
downY = max(downY, 0);
        upY = max(upY, 0);
    }

    int2 cellCoord0 = int2(leftX, downY);
    int2 cellCoord1 = int2(leftX, upY);
    int2 cellCoord2 = int2(rightX, upY);
    int2 cellCoord3 = int2(rightX, downY);

    CellData cellData0 = _CellDatas[CellCoordToIndex(cellCoord0)];
    CellData cellData1 = _CellDatas[CellCoordToIndex(cellCoord1)];
    CellData cellData2 = _CellDatas[CellCoordToIndex(cellCoord2)];
    CellData cellData3 = _CellDatas[CellCoordToIndex(cellCoord3)];

    float tempDensity0 = lerp(cellData0.density, cellData1.density, vt);
    float tempDensity1 = lerp(cellData3.density, cellData2.density, vt);
    float finalDensity = lerp(tempDensity0, tempDensity1, ut);

    float2 tempVelocity0 = lerp(cellData0.velocity, cellData1.velocity, vt);
    float2 tempVelocity1 = lerp(cellData3.velocity, cellData2.velocity, vt);
    float2 finalVelocity = lerp(tempVelocity0, tempVelocity1, ut);

    float4 tempColor0 = lerp(cellData0.color, cellData1.color, vt);
    float4 tempColor1 = lerp(cellData3.color, cellData2.color, vt);
    float4 finalColor = lerp(tempColor0, tempColor1, ut);

    AdvectData advectData = _AdvectDatas[index];
    advectData.density = finalDensity;
    advectData.velocity = finalVelocity;
    advectData.color = finalColor;
    _AdvectDatas[index] = advectData;
}

#pragma kernel CSMain

#include "Assets/FluidSimulationLibrary.hlsl"

RWStructuredBuffer<CellData> _CellDatas;
RWStructuredBuffer<AdvectData> _AdvectDatas;

[numthreads(256,1,1)]
void CSMain (uint3 id : SV_DispatchThreadID)
{
    uint index = id.x;
    if(index >= _CellDatas.Length) {
        return;
    }

    CellData cellData = _CellDatas[index];
    AdvectData advectData = _AdvectDatas[index];

    cellData.density = advectData.density;
    cellData.velocity = advectData.velocity;
    cellData.color = advectData.color;

    advectData.density = 0;
    advectData.velocity = 0;
    advectData.color = 0;

    _CellDatas[index] = cellData;
    _AdvectDatas[index] = advectData;
}

衰减

#pragma kernel CSMain

#include "Assets/FluidSimulationLibrary.hlsl"

RWStructuredBuffer<CellData> _CellDatas;

float _DensityDamping;
float _VelocityDamping;

[numthreads(256,1,1)]
void CSMain (uint3 id : SV_DispatchThreadID)
{
    uint index = id.x;
    if(index >= _CellDatas.Length) {
        return;
    }

    CellData cellData = _CellDatas[index];
    cellData.density *= _DensityDamping;
    cellData.velocity *= _VelocityDamping;
    cellData.color *= _DensityDamping;
    _CellDatas[index] = cellData;
}

现在，已经有了一个基础效果：

涡度约束（Vorticity Confinement）
涡度约束的作用是向流体加入卷曲的运动趋势，让流体运动更符合自然规律。

struct VortexData {
    float2 velocity;
};

#pragma kernel CSMain

#include "Assets/FluidSimulationLibrary.hlsl"

StructuredBuffer<CellData> _CellDatas;
RWStructuredBuffer<VortexData> _VortexDatas;

float _VortexIntensity;

float GetCurl(int2 coord) {
    int2 leftCellCoord = coord + int2(-1, 0);
    int2 rightCellCoord = coord + int2(1, 0);
    int2 upCellCoord = coord + int2(0, 1);
    int2 downCellCoord = coord + int2(0, -1);

    int leftCellIndex = CellCoordToIndex(leftCellCoord);
    int rightCellIndex = CellCoordToIndex(rightCellCoord);
    int upCellIndex = CellCoordToIndex(upCellCoord);
    int downCellIndex = CellCoordToIndex(downCellCoord);

    CellData leftCellData = _CellDatas[leftCellIndex];
    CellData rightCellData = _CellDatas[rightCellIndex];
    CellData upCellData = _CellDatas[upCellIndex];
    CellData downCellData = _CellDatas[downCellIndex];

    return upCellData.velocity.x - downCellData.velocity.x + leftCellData.velocity.y - rightCellData.velocity.y;
}

[numthreads(256,1,1)]
void CSMain (uint3 id : SV_DispatchThreadID)
{
    uint index = id.x;
    if(index >= _CellDatas.Length) {
        return;
    }

    CellData cellData = _CellDatas[index];
    if(cellData.coord.x < 2 || cellData.coord.x > _Resolution.x - 3 || cellData.coord.y < 2 || cellData.coord.y > _Resolution.y - 3) {
        return;
    }

    int2 leftCellCoord = cellData.coord + int2(-1, 0);
    int2 rightCellCoord = cellData.coord + int2(1, 0);
    int2 upCellCoord = cellData.coord + int2(0, 1);
    int2 downCellCoord = cellData.coord + int2(0, -1);

    float centerCurl = GetCurl(cellData.coord);
    float leftCurl = GetCurl(leftCellCoord);
    float rightCurl = GetCurl(rightCellCoord);
    float upCurl = GetCurl(upCellCoord);
    float downCurl = GetCurl(downCellCoord);

    float dx = abs(downCurl) - abs(upCurl);
    float dy = abs(rightCurl) - abs(leftCurl);
    float len = sqrt(dx * dx + dy * dy);
    if(len == 0) {
        return;
    }

    dx = _VortexIntensity / len * dx;
    dy = _VortexIntensity / len * dy;

    float scaler = length(cellData.velocity) * saturate(cellData.density * 10);

    VortexData vortexData = _VortexDatas[index];
    vortexData.velocity.x += centerCurl * dx * scaler;
    vortexData.velocity.y += centerCurl * dy * scaler;
    _VortexDatas[index] = vortexData;
}

#pragma kernel CSMain

#include "Assets/FluidSimulationLibrary.hlsl"

RWStructuredBuffer<CellData> _CellDatas;
RWStructuredBuffer<VortexData> _VortexDatas;

float _MaxSpeed;

[numthreads(256,1,1)]
void CSMain (uint3 id : SV_DispatchThreadID)
{
    uint index = id.x;
    if(index >= _VortexDatas.Length) {
        return;
    }

    CellData cellData = _CellDatas[index];
    VortexData vortexData = _VortexDatas[index];

    cellData.velocity += vortexData.velocity;

    float speed = length(cellData.velocity);
    if(speed > _MaxSpeed) {
        cellData.velocity = normalize(cellData.velocity) * _MaxSpeed;
    }

    vortexData.velocity = 0;

    _CellDatas[index] = cellData;
    _VortexDatas[index] = vortexData;
}

加入涡度约束后，效果更自然了：

以上即为流体模拟的主要计算过程。

五、结语

得益于显卡的快速发展，已经有PC游戏开始使用实时流体模拟了。相对于传统的粒子特效，用流体模拟做烟、云这类效果，最大的优势就是可交互性强，角色从烟雾中穿过，烟雾会被拨开，飞机从云层中穿过，云会被冲散。Houdini里惊艳的影视特效，很多就是用流体模拟的方法实现的。

这次研究流体模拟的初衷，是想尝试在Unity里做一个流体特效引擎，现在只完成了最基础的2D模拟，距离最终目标还很遥远。

源文件
Github：https://github.com/MagicStones23/Unity-Fluid-Simulation

百度网盘：https://pan.baidu.com/s/14kqkyxjikb3cguN55y_X7w?pwd=1111
提取码：1111

这是侑虎科技第1507篇文章，感谢作者异世界的魔法石供稿。欢迎转发分享，未经作者授权请勿转载。如果您有任何独到的见解或者发现也欢迎联系我们，一起探讨。（QQ群：465082844）

作者主页：https://www.zhihu.com/people/shui-guai-76-84

再次感谢异世界的魔法石的分享，如果您有任何独到的见解或者发现也欢迎联系我们，一起探讨。（QQ群：465082844）

暂无标签