In this article, you will learn about fluid simulation using warping, backward advection and assuming a periodic velocity.
Introduction
This article takes a look at Jos Stams Stable Fluids, the current standard for fluid simulation in games and beyond, and proposes a less complicated but also less scientific oriented approach.
Background
What I was looking for was a paint stirring simulation. I came across Jos Stam's Stable Fluids. Its implementation does look like stirred fluids with high viscosity. It is a very popular method for fluid animation and is supposed to be scientifically sound as well, solving the Navier-Stokes equations. Alas, I do not understand it. Originally, it was presented as a method for fast fluid simulation for games, but many scientists make mention of it now. The first code was about injecting smoke into a liquid tank, the liquid an incompressable fluid. And I have three questions here.
- If I inject smoke into a liquid, it will bubble up, not spread like when I stir 2 colors of paint.
- If the tank is filled to the brim with an incompressible liquid, adding volume to it will make the tank flow over or explode.
- The implementation of advection. Advection generally means transport of mass, heat, momentum. The specific form here is a backward tracking.
What follows here is my simple, basic, experiments with fluid flow simulation and backward tracking advection. I do not have much knowledge of fluid mechanics. My implementation has no scientific value whatsoever. It is merely an effort to simulate fluid flow.
This is my simulation of the Kelvin-Helmholtz instability:
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
Vel(tX, tY, 2) = qy1
For tX2 = 0 To 2
ImageData(tX2, tX, tY) = 255
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
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
Intfield(tX, tY, 2) = Intfield(tX1, tY1, 1)
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
It starts with a white field above a black, and I assign velocities. Not a surprise that the top velocity is opposite to the bottom one. But I also assume a periodical, sinusoid velocity perpendicular (y-direction) to these x-direction velocities. That is not supported by science. It was just something that I noticed in many pictures and videos of real fluid movement: Rayleigh-Taylor instability (looks like a variation to Kelvin-HelmHoltz), ink in water (looks like a form of Rayleigh-Taylor), and incense smoke columns. Some of these do not look like the Stable Fluid simulations. It looks as if there is a periodic velocity as well. A bit like low frequency sound waves, with a velocity straight ahead in horizontal direction accompanied by periodic amplitudes in vertical direction. It seems like something I can use in my fluid simulation, without explaining how it was caused, just assuming it.
I originally used this line to assign a value to every pixel:
Intfield(tX, tY, 2) = Intfield(tX1, tY1, 1)
And the result is a curly wave. If I use the straightforward (literally) variation:
Intfield(tX1, tY1, 2) = Intfield(tX, tY, 1)
I get no curls! The curl in the picture is caused by the backward advection.
What happens is that a black pixel occasionally taps into the velocity of the white field, and vice versa.
It did not really come from there but the backward calculation makes it look that way. Backward advection causes all those nice curls. (I have seen many implementations where extra curl or vorticity is added. That does make even nicer pictures.)
As I said, I was originally looking for a simulation of stirring 2 colors of paint. It can be simulated by warping. That is the Liquify tool in Photoshop. Here is my implementation:
Private Sub InitWarp
ReDim cSeg1(-rS To rS, -rS To rS) As Single
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))
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
dy = yM - y
qX5 = pBox1.Width
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(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)
Next t
End If
End If
Next qX
End If
Next qY
End Sub
Then when I looked at soap bubbles, the movements of the colored soap film, it seemed like an all-over warp. I implemented it like this:
Private Sub InitRand(b1 As Integer, b2 As Integer)
Randomize
Dim Randi(0 To 511) As Integer
tY1 = 255
tX1 = CInt(Rnd * 2)
If tX1 < 1 Then tX1 = -1
Randi(0) = tX1
For tX = 1 To 511
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
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
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
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
The drawing routine is the same as above. Again, the backward advection creates curls and fingers.
Further Reading
Conclusion and Points of Interest
Decide for yourself if this very simple algorithm creates realistic fluid simulations! It was coded in VB6, but by using C++ and a GPU implementation it can be much faster. I hope to inspire you to experiment more. As for myself, next I want to try to get better simulations of smoke columns and ink in water, as well as soap bubble films. Suggestions are welcome!
History
- 10th November, 2020: Initial version
My name is Franciska Ruessink. I studied IT at The Hague Hogeschool (University) and was employed in IT as software developer, database developer, network administrator, and IT project manager for several years.
After an early retirement I continue to write software for the fun of it and to help friends, I also did some web development. I usually work with VB but started again with VC++ recently.