【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版本是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://pan.baidu.com/s/14kqkyxjikb3cguN55y_X7w?pwd=1111
提取码：1111
这是侑虎科技第1507篇文章，感谢作者异世界的魔法石供稿。欢迎转发分享，未经作者授权请勿转载。如果您有任何独到的见解或者发现也欢迎联系我们，一起探讨。（QQ群：465082844）
作者主页：
再次感谢异世界的魔法石的分享，如果您有任何独到的见解或者发现也欢迎联系我们，一起探讨。（QQ群：465082844）
,
