首页 > 其他分享 >技术美术｜游戏中的流体模拟（Fluid Simulation）

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

时间：2023-12-15 15:22:36浏览次数：31  
标签：index cellData Fluid staggeredPoint Simulation 流体 pointCoord int2 velocity

【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）

标签：index,cellData,Fluid,staggeredPoint,Simulation,流体,pointCoord,int2,velocity
From： https://www.cnblogs.com/uwatech/p/17903440.html

相关文章

  • How to Use Docker and NS-3 to Create Realistic Network Simulations
    https://insights.sei.cmu.edu/blog/how-to-use-docker-and-ns-3-to-create-realistic-network-simulations/ HowtoUseDockerandNS-3toCreateRealisticNetworkSimulationsALEJANDROGOMEZMARCH27,2023Sometimes,researchersanddevelopersneedt......
  • 多开器在Windows电脑上的流体力学仿真应用
    多开器在Windows电脑上的流体力学仿真应用摘要：随着计算机技术的不断发展，流体力学仿真成为了研究和解决涉及液体和气体运动问题的重要工具。而在Windows电脑上，多开器的出现为流体力学仿真应用提供了更大的便利性和效率。本文将介绍多开器在Windows电脑上的流体力学仿真应用，并探讨......
  • 在fluid主题中加入Google广告
    title:在fluid主题中加入Google广告banner_img:https://proxy.thisis.plus/8592ed575a242368611755f5529c28e.pngdate:2023-1-2710:00:00categories:-踩坑在fluid主题中加入Google广告在fluid的官方文档中，提供了在fluid主题中加入Google广告的方法，但是其中提到的参数d......
  • DoraOS云终端moonlight串流体验
    前言大家如果手里有这种需求，主机在书房，想在客厅或者卧室使用主机的强大性能，搬动主机又很麻烦。这时，串流就可以完美的解决这个问题，随便一个老旧pc或者瘦客户机，安装DoraOS软件，即可实现低延时，高画质的串流，在客厅看电影，插上手柄玩主机游戏。串流是一种数字传输方式，它可以将大......
  • 首个流体力学大模型背后，是昇腾的大模型“造林”逻辑
    作者|曾响铃文|响铃说一个飞机模型在试验风洞里，空气从它的机翼与机身流过，形成一层又一层稳定的气流，当风速加快，空气的流线开始波浪式摆动，最终随着速度增大而相互混合、形成不再能分辨的湍流，看起来混沌又无序……这是流体力学测试的常见场景，一遍又一遍地测试，只为模拟或预测真实的......
  • 应用动量定理处理流体问题
    建立流体模型对于一段流体质量具有连续性，其密度为\(ρ\)流速为\(v\)流体横截面积为\(S\)微元研究微元作用时间：\(Δt\)微元作用长度：\(vΔt\)则对应的质量为：\[Δm=ρSvΔt\]随后建立方程，应用动量定理研究即可。......
  • uvm获取simulation option
    获取simoption的传统方式是用$test$plusargs(),和$value$plusargs().UVM提供了uvm_cmdline_processor来获取simopt.用法如下：uvm_cmdline_processorclp=uvm_cmdline_processor::get_inst();stringval_q[$];clp.get_arg_values("+val=",val_q);uvm_cmdline_processor......
  • Solidworks流体仿真插件安装及案例分析
    Solidworks流体仿真插件安装及案例分析1流体仿真插件的安装如图1所示，安装时勾选SolidworksFlowSimuation模块，一路“下一步”安装完毕。完成安装后打开软件，图2所示，点击Solidworks插件按钮，找到SolidworksFlowSimuation按钮，双击可打开说明插件安装成2案例分析2.1案例背景......
  • 基于 ACK Fluid 的混合云优化数据访问（四）：将第三方存储目录挂载到 Kubernetes，提升效率和
    作者：车漾前文回顾：本系列将介绍如何基于ACKFluid支持和优化混合云的数据访问场景，相关文章请参考：-基于ACKFluid的混合云优化数据访问（一）：场景与架构-基于ACKFluid的混合云优化数据访问（二）：搭建弹性计算实例与第三方存储的桥梁-基于ACKFluid的混合云优化数据访问（三）：加速......
  • 基于 ACK Fluid 的混合云优化数据访问（三）：加速第三方存储的读访问，降本增效并行
    作者：车漾前文回顾：本系列将介绍如何基于ACKFluid支持和优化混合云的数据访问场景，相关文章请参考：基于ACKFluid的混合云优化数据访问（一）：场景与架构基于ACKFluid的混合云优化数据访问（二）：搭建弹性计算实例与第三方存储的桥梁在前一篇文章《搭建弹性计算实例与第三方存储的......

赞助商

阅读排行

网站主页关于我们联系我们网站地图

本网站内容转载自其他媒体，侵权联系[admin##ips99.com]。

Copyright © 2020-2023 IPS99 版权所有 IPS99