On Fluid Simulation

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.

1. If I inject smoke into a liquid, it will bubble up, not spread like when I stir 2 colors of paint.
2. If the tank is filled to the brim with an incompressible liquid, adding volume to it will make the tank flow over or explode.
3. 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 ''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

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) ''backward tracking advection

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 ''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

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 ''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

The drawing routine is the same as above. Again, the backward advection creates curls and fingers.

Further Reading

• Presentation by Jos Stam https://www.ljll.math.upmc.fr/~frey/cours/references/Stam%20J.,%20Stable%20fluids.pdf
  Many more on the subject, including source codes for various languages if you type "Jos Stam Stable Fluids" in your search engine.
• Kevin-Helmholtz instability: https://en.wikipedia.org/wiki/Kelvin%E2%80%93Helmholtz_instability
• Rayleigh-Taylor instability: https://en.wikipedia.org/wiki/Rayleigh%E2%80%93Taylor_instability
• Ink in water: https://www.youtube.com/watch?v=t-Are2N7zUE
• Incense smoke: https://photodune.net/item/abstract-shapes-formed-by-scrolls-of-incense-smoke/21529834 (very visible are the periodic rings that start at the bottom of the smoke column as increase and decrease of the smoke density, and become curls higher up the smoke column)

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

• 10^th November, 2020: Initial version