16,000,304 members
Articles / Programming Languages / Visual Basic 6

# On Fluid Simulation

Rate me:
12 Nov 2020CPOL4 min read 6.4K   4   1
A new approach to 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:

VB.NET
```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:

VB.NET
`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:

VB.NET
`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:

VB.NET
```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:

VB.NET
```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.

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

Written By
Retired
Netherlands
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.