技术美术|游戏中的流体模拟(Fluid Simulation)
【USparkle专栏】如果你深怀绝技,爱“搞点研究”,乐于分享也博采众长,我们期待你的加入,让智慧的火花碰撞交织,让知识的传递生生不息!
一、闲聊
最近一直在研究流体模拟,很神奇的一个东西,好在网上有很多参考资料,研究过程不算太困难。分享下最近一段时间的学习心得。
二、效果演示
三、算法原理
游戏领域实现流体模拟的几种常见方式有:
- 基于网格的方法:在网格上模拟,每个格子都有自己的数据(速度、密度、颜色、温度等),逐帧更新格子内数据。这种方法的优点是方便多线程实现,渲染也很方便。缺点是计算过程中需要对参数做估算,容易产生误差。
- 基于粒子的方法:将流体具象化为很多个粒子,每个粒子都有自己数据(速度、颜色、温度),逐帧更新粒子的位置。这种方法的优点是误差小,能表现出更多的流体细节。缺点是不利于多线程实现,渲染也比较麻烦。
这篇文章采用的是基于网格的方法,流体有很多种类(气体、水、岩浆、蜂蜜等),不同流体使用的算法各有差异,这篇文章讨论的是气体流体模拟。
在流体模拟中,有两个主要计算过程,压缩解算和流动。
压缩解算(Project)
压缩性是流体的基本属性之一,正常环境下,大多数流体都很难被压缩,向流体施加很大的力,而流体的体积变化却很小。
压缩解算的目的,就是要模拟流体很难被压缩的特点,假设我们在一个8x8的网格上做流体模拟:

先在格子边框上创建辅助点(Staggered Point),水平方向辅助点为黄色,垂直方向辅助点为绿色:

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

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

单个格子拆分前

单个格子拆分后

整体拆分前

整体拆分后
然后根据格子周围4个辅助点的速度,对格子做压缩解算:

上图的格子有三个方向在流入,一个方向在流出,流入量大于流出量,要使流体不被压缩,流入量和流出量必须相等。
先计算净流入、流出量(Divergence):
1 | float divergence = leftPointSpeed - rightPointSpeed + downPointSpeed - upPointSpeed; |
将其均分后修改辅助点速度:
1 2 3 4 5 | divergence /= 4; leftPointSpeed += -divergence; rightPointSpeed += divergence; downPointSpeed += -divergence; upPointSpeed += divergence; |

修改辅助点速度后
这样就保证了这个格子流入和流出量相等。
再看一个流体遇到障碍物的例子:

格子的右边是一面墙,所以右边黄色点的速度始终为0,压缩解算的公式变为:
1 2 3 4 5 6 | float divergence = leftPointSpeed + downPointSpeed - upPointSpeed; divergence /= 3; leftPointSpeed += -divergence; downPointSpeed += -divergence; upPointSpeed += divergence; |

修改辅助点速度后
我们可以为每个辅助点附加一个Scaler,障碍物的Scaler为0,非障碍物的Scaler为1,这样一来,有无障碍物都可以使用同一个公式:
1 2 3 4 5 6 7 8 9 10 11 12 13 | 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; |
最后,计算辅助点速度的平均值,更新格子速度:
1 2 3 4 | float uSpeed = (leftPointSpeed + rightPointSpeed) / 2; float vSpeed = (downPointSpeed + upPointSpeed) / 2; cellData.speed = float2(uSpeed, vSpeed); |
一个格子速度发生变化,其临近格子的流入、流出量也会改变,这里我们需要迭代多次去逼近正确解:
1 2 3 4 5 6 7 | for ( int i = 0; i < iteration; i++) { for ( int j = 0; j < cellNumber; j++) { //对格子做计算,修改辅助点速度 } //更新格子速度 } |
这样,压缩解算就完成了!
流动(Advect)
更新完格子的速度后,就可以移动格子内的数据了。最直观的做法是,根据格子的速度,计算出它移动到了哪个位置,然后把它的数据(密度,速度等)加入到新格子中。

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

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

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

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

扩散(Diffuse)
流体还有扩散特性,高浓度区域会主动扩散到低浓度区域,直至所有格子的浓度相等。比方说向水杯里滴入一滴墨水,墨水会逐渐扩散开,直至整杯水均匀变黑。不过这一步并不是必需的,我在实际尝试中发现,加入扩散后效果反而没那么好看了。
四、Unity内具体实现过程
我使用的Unity版本是2021.3,URP管线。流体模拟的计算量比较大,我使用的ComputeShader做计算。
主要流程
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 | private void OnEnable() { //定义并初始化数据结构 } private void Update() { //向指定格子输入流体 //将格子速度拆分到辅助点上 //压缩解算,修改辅助点速度 //更新格子速度 //格子数据流动 //衰减 } |
数据结构
CellData为单个格子的数据结构,UStaggeredPoint代表水平方向的辅助点,VStaggeredPoint代表垂直方向的辅助点。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 | 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 |
初始化
初始化格子数据:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 | #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两次,分别初始化水平和垂直方向的辅助点:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 | #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,找到鼠标周围的格子,修改数据:
1 2 3 4 5 6 7 | struct InjectData { int2 center; float radius; float density; float2 velocity; float4 color; }; |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 | #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; } |
速度拆分
将格子的速度拆分到辅助点上,先累加,再求平均:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 | #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; } } |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 | #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; } } |
压缩解算
根据净流入、流出量修改辅助点速度,先累加,再求平均:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 | #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; } |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 | #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; } } |
更新格子速度
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 | #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个格子内的数据做插值,更新自己:
1 2 3 4 5 | struct AdvectData { float density; float2 velocity; float4 color; }; |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 | #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; } |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 | #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; } |
衰减
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 | #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)
涡度约束的作用是向流体加入卷曲的运动趋势,让流体运动更符合自然规律。
1 2 3 | struct VortexData { float2 velocity; }; |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 | #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; } |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 | #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)
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 清华大学推出第四讲使用 DeepSeek + DeepResearch 让科研像聊天一样简单!
· 推荐几款开源且免费的 .NET MAUI 组件库
· 实操Deepseek接入个人知识库
· 易语言 —— 开山篇
· Trae初体验
2022-12-15 关于切换场景加载耗时的优化问题
2021-12-15 在UI上制作动画的方案选择
2020-12-15 带BlendShape表情的动作文件播放异常
2020-12-15 如何通过Timeline的形式实现技能编辑器