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

大家好，欢迎来到IT知识分享网。

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

一、闲聊

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

二、效果演示

三、算法原理

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

• 基于网格的方法：在网格上模拟，每个格子都有自己的数据（速度、密度、颜色、温度等），逐帧更新格子内数据。这种方法的优点是方便多线程实现，渲染也很方便。缺点是计算过程中需要对参数做估算，容易产生误差。

• 基于粒子的方法：将流体具象化为很多个粒子，每个粒子都有自己数据（速度、颜色、温度），逐帧更新粒子的位置。这种方法的优点是误差小，能表现出更多的流体细节。缺点是不利于多线程实现，渲染也比较麻烦。

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

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

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

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

先在格子边框上创建辅助点（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：GitHub – MagicStones23/Unity-Fluid-Simulation

百度网盘：百度网盘 请输入提取码
提取码：1111

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

作者主页：异世界的魔法石 – 知乎

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