关于流体模拟
5.00/5 (3投票s)
一种新的流体模拟方法
引言
本文探讨了 Jos Stam 的稳定流体，这是当前游戏及其他领域流体模拟的标准，并提出了一种不那么复杂但科学性也较低的方法。
背景
我最初想寻找的是一个搅拌颜料的模拟。我遇到了 Jos Stam 的稳定流体。它的实现看起来确实像是高粘度的搅拌流体。这是一种非常流行的流体动画方法，并且据称在科学上也是站得住脚的，它解决了 Navier-Stokes 方程。唉，我不明白它。最初，它被提出作为游戏快速流体模拟的一种方法，但现在许多科学家都提到了它。第一个代码是关于将烟雾注入液体罐中，该液体是不可压缩的流体。在这里我有三个问题。
- 如果我将烟雾注入液体中，它会冒泡上升，而不是像我搅拌两种颜色的颜料那样扩散。
- 如果罐子装满了不可压缩的液体，向其中添加体积会使罐子溢出或爆炸。
- 对流的实现。对流通常意味着质量、热量、动量的输送。这里的特定形式是后向追踪。
接下来是我对流体流动模拟和后向追踪对流的简单、基本的实验。我对流体力学知之甚少。我的实现没有任何科学价值。它仅仅是为了模拟流体流动。
这是我对 Kelvin-Helmholtz 不稳定性的模拟
Private Sub InitiateHelm()
pi16 = Pi / 48
For tX = 0 To 511
qy1 = Sin(tX * pi16) * tX * 0.01
For tY = 0 To 255
Vel(tX, tY, 1) = 1 ''velocity x
Vel(tX, tY, 2) = qy1 ''velocity y
For tX2 = 0 To 2
ImageData(tX2, tX, tY) = 255 ''white top half of the picture
Next tX2
Intfield(tX, tY, 1) = 255
Next tY
For tY = 256 To 511
Vel(tX, tY, 1) = -1
Vel(tX, tY, 2) = qy1
For tX2 = 0 To 2
ImageData(tX2, tX, tY) = 0 ''black bottom half of the picture
Next tX2
Intfield(tX, tY, 1) = 0
Next tY
Next tX
Private Sub Draw()
Dim a1 As Double, b1 As Double, qx1 As Double, qy1 As Double
For tX = 0 To 510
For tY = 0 To 510
qx1 = tX + Vel(tX, tY, 1)
qy1 = tY + Vel(tX, tY, 2)
tX1 = Int(qx1)
tY1 = Int(qy1)
qx1 = qx1 - tX1
qy1 = qy1 - tY1
If tX1 < 0 Then tX1 = 0
If tX1 > 510 Then tX1 = 510
If tY1 < 0 Then tY1 = 0
If tY1 > 510 Then tY1 = 510
''a1 = lerp(CDbl(Intfield(tX1, tY1, 1)),
''CDbl(Intfield(tX1 + 1, tY1, 1)), qx1) ''interpolation for smoother picture
''b1 = lerp(CDbl(Intfield(tX1, tY1 + 1, 1)), CDbl(Intfield(tX1 + 1, tY1 + 1, 1)), qx1)
''Intfield(tX, tY, 2) = lerp(a1, b1, qy1)
Intfield(tX, tY, 2) = Intfield(tX1, tY1, 1) ''backward tracking advection
''Intfield(tX1, tY1, 2) = Intfield(tX, tY, 1) ''plain forward
Next tY
Next tX
For tX = 0 To 511
For tY = 0 To 511
Intfield(tX, tY, 1) = (Intfield(tX , tY, 2)
For tX2 = 0 to 2
ImageData(tX2, tX, tY) = Intfield(tX, tY, 1)
Next tx2
Next tY
Next tX
End Sub
Using the Code
它从一个白色的场上方开始，我分配了速度。顶部的速度与底部的速度相反，这并不奇怪。但我也假设一个周期性的、正弦的速度，该速度垂直（y 方向）于这些 x 方向的速度。这在科学上是不成立的。这只是我在许多真实流体运动的图片和视频中注意到的现象：Rayleigh-Taylor 不稳定性（看起来像是 Kelvin-HelmHoltz 的一种变体），墨水在水中（看起来像是 Rayleigh-Taylor 的一种形式），以及香烟烟柱。其中一些看起来不像稳定流体模拟。它看起来好像也有一个周期性的速度。有点像低频声波，在水平方向上有一个直线前进的速度，伴随着垂直方向上的周期性振幅。这似乎是我可以在我的流体模拟中使用的东西，而无需解释它是如何造成的，只是假设它。
我最初使用此行来为每个像素赋值
Intfield(tX, tY, 2) = Intfield(tX1, tY1, 1) ''backward tracking advection
结果是一个卷曲的波浪。如果我使用直接（字面意义上）的变体
Intfield(tX1, tY1, 2) = Intfield(tX, tY, 1)
我没有得到任何卷曲！图片中的卷曲是由后向对流引起的。
发生的事情是，一个黑色的像素偶尔会利用白色场的速度，反之亦然。
它实际上不是来自那里，但后向计算让它看起来是这样的。后向对流导致了所有那些漂亮的卷曲。（我见过许多实现，其中添加了额外的卷曲或涡度。这确实能做出更好的图片。）
正如我所说，我最初想寻找的是搅拌两种颜色的颜料的模拟。这可以通过变形来模拟。这就是 Photoshop 中的液化工具。这是我的实现
Private Sub InitWarp
ReDim cSeg1(-rS To rS, -rS To rS) As Single ''rS is the radius of the warp brush
For tX = 0 To rS
dax = tX * tX
For tY = 0 To rS
dxy = Sqr(dax + tY * tY)
If dxy <= rS Then
dr = LineArr(Int(dxy)) ''LineArr stores the warp distances in the brush,
''larger distance of warp = closer to the middle of the brush
cSeg1(tX, tY) = dr
cSeg1(-tX, tY) = dr
cSeg1(tX, -tY) = dr
cSeg1(-tX, -tY) = dr
End If
Next tY
Next tX
End Sub
Private Sub DoWarp
dx = xM - x ''xM = middle(old mouse position x), x = new mouse position x
dy = yM - y
qX5 = pBox1.Width ''Actual picturebox
qY5 = pBox1.Height
For qY = -rS To rS
qY1 = qY + yM
If qY1 > -1 And qY1 < qY5 Then
For qX = -rS To rS
qX1 = qX + xM
If qX1 > -1 And qX1 < qX5 Then
ShiftArr(qX1, qY1, 1) = ShiftArr(qX1, qY1, 1) + cSeg1(qX, qY) * dx ''Shiftarr is a
''buffer that stores the warp distances of the brush, applied to the picture
ShiftArr(qX1, qY1, 2) = ShiftArr(qX1, qY1, 2) + cSeg1(qX, qY) * dy
End If
Next qX
End If
Next qY
For qY = -rS To rS
qY1 = qY + yM
If qY1 > -1 And qY1 < qY5 Then
For qX = -rS To rS
qX1 = qX + xM
If qX1 > -1 And qX1 < qX5 Then
qX3 = qX1 + ShiftArr(qX1, qY1, 1)
qY3 = qY1 + ShiftArr(qX1, qY1, 2)
If qX3 > -1 And qX3 < qX5 And qY3 > -1 And qY3 < qY5 Then
For t = 0 To 2
ImageData(t, qX1, qY1) = ImageData2(t, qX3, qY3) ''Imagedata2 is the picture,
''Imagedata is a copy that will replace Imagedata2 once all changes are drawn
Next t
End If
End If
Next qX
End If
Next qY
End Sub
然后，当我观察肥皂泡时，彩色肥皂膜的运动，它看起来像一个整体的变形。我像这样实现了它
Private Sub InitRand(b1 As Integer, b2 As Integer)
Randomize
Dim Randi(0 To 511) As Integer ''random warp distance across the picture
tY1 = 255
tX1 = CInt(Rnd * 2)
If tX1 < 1 Then tX1 = -1
Randi(0) = tX1
For tX = 1 To 511 ''y direction, random decide to add or distract 1
''from the previous warp distance
qy1 = Rnd
If qy1 < 0.167 Then tX1 = Abs(tX1) - 1
If qy1 > 0.833 Then tX1 = -Abs(tX1) + 1
If Randi(tX - 1) < -b1 Then tX1 = 1
If Randi(tX - 1) > b1 Then tX1 = -1
Randi(tX) = Randi(tX - 1) + tX1
Next tX
For tX = 3 To 508 ''smoothen the random values a bit
Vel(tX, 256, 2) = (Randi(tX) + Randi(tX - 1) + Randi(tX + 1) + _
Randi(tX - 2) + Randi(tX + 2) + Randi(tX - 3) + Randi(tX + 3)) / 7
Next tX
For tX = 0 To 511
qx1 = Vel(tX, 256, 2) / 256
For tY = 0 To 255
Vel(tX, tY, 2) = tY * qx1 ''warp distance at the edges = 0,
'' calculate the warp (velocity) distance for all pixels, y direction
Next tY
qx1 = -qx1
For tY = 1 To 255
Vel(tX, tY + 256, 2) = Vel(tX, 256, 2) + tY * qx1
Next tY
qy1 = Vel(tX, 511, 2)
Next tX
tY1 = 255
tX1 = CInt(Rnd * 2)
If tX1 < 1 Then tX1 = -1
Randi(0) = tX1
For tX = 1 To 511 '' same for x direction
qy1 = Rnd
If qy1 < 0.167 Then tX1 = Abs(tX1) - 1
If qy1 > 0.833 Then tX1 = -Abs(tX1) + 1
If Randi(tX - 1) < -b2 Then tX1 = 1
If Randi(tX - 1) > b2 Then tX1 = -1
Randi(tX) = Randi(tX - 1) + tX1
Next tX
For tX = 3 To 508
Vel(256, tX, 1) = (Randi(tX) + Randi(tX - 1) + Randi(tX + 1) + _
Randi(tX - 2) + Randi(tX + 2) + Randi(tX - 3) + Randi(tX + 3)) / 7
Next tX
For tX = 0 To 511
qx1 = Vel(256, tX, 1) / 256
For tY = 0 To 255
Vel(tY, tX, 1) = tY * qx1
Next tY
qx1 = -qx1
For tY = 1 To 255
Vel(tY + 256, tX, 1) = Vel(256, tX, 1) + tY * qx1
Next tY
Next tX
End Sub
绘图程序与上面相同。同样，后向对流产生卷曲和手指。
延伸阅读
- Jos Stam 的演示 https://www.ljll.math.upmc.fr/~frey/cours/references/Stam%20J.,%20Stable%20fluids.pdf
更多关于这个主题的内容，包括各种语言的源代码，如果您在搜索引擎中键入“Jos Stam Stable Fluids”。 - Kevin-Helmholtz 不稳定性：https://en.wikipedia.org/wiki/Kelvin%E2%80%93Helmholtz_instability
- Rayleigh-Taylor 不稳定性：https://en.wikipedia.org/wiki/Rayleigh%E2%80%93Taylor_instability
- 墨水在水中：https://www.youtube.com/watch?v=t-Are2N7zUE
- 香烟烟雾：https://photodune.net/item/abstract-shapes-formed-by-scrolls-of-incense-smoke/21529834（非常明显的是从烟柱底部开始的周期性环，烟雾密度增加和减少，并在烟柱上方变成卷曲）
结论和关注点
自己决定这种非常简单的算法是否能创建逼真的流体模拟！它使用 VB6 编码，但通过使用 C++ 和 GPU 实现，它可以更快。我希望激励您进行更多实验。至于我自己，接下来我想尝试更好地模拟烟柱和墨水在水中的情况，以及肥皂泡膜。欢迎提出建议！
历史
- 2020 年 11 月 10 日：初始版本