In this tip, you will learn about the elements of many current fluid simulation algorithms that make the fluid velocity curl.
Introduction
Fluid simulation is fascinating. Fluid moves in wonderful curves, we can see it when mixing different colors of fluid. I came by the subject because of my love for Art Nouveau, or Liberty style. There the lines in designs seem to be inspired by fluid motions. There are a lot of fluid simulation algorithms. Fluid has velocity, in a 2D simulation, it is a vector, consisting of quantities along the x-axis and y-axis. That means 2 arrays, u and v, that hold the velocity in x-direction and y direction for every point in the plane. But fluid does not just move in a straight line. It ondulates. Let's take a look at how the algorithms make the velocity curl.
Projection
Code Step 1
Curling velocity is achieved by adding some of the x-direction velocity to the y-direction, and vice versa. The idea behind it is that the velocity at every point is accompanied by a secondary velocity perpendicular to it. So part of the velocity is rotated at an angle, the result is that the movement of the fluid is not in a straight line, but it curls. Algorithms usually calculate this with the help of an intermediate array, some call it pressure, others vorticity. The most basic manipulation (in a N1 x N1 plane) is merging u and v into pressure array p and then splitting it again:
k = -0.5
For i = 1 to 10
For tX = 0 To N1
For tY = 0 To N1
p(tX, tY) = k * (u(tX + 1, tY) - u(tX - 1, tY) + (v(tX, tY + 1) - v(tX, tY - 1)))
Next tY
Next tX
For tX = 1 To N
For tY = 1 To N
u(tX, tY) = u(tX, tY) + k * (p(tX + 1, tY) - p(tX - 1, tY))
v(tX, tY) = v(tX, tY) + k * (p(tX, tY + 1) - p(tX, tY - 1))
Next tY
Next tX
next i
This produces curly velocities.
Code Step 2
Taking it one step further, the intermediate array can be manipulated between the merging and splitting step to produce extra curl. This is done by applying linear solving. The right-hand side rhs gets the merged values from u and v, that is the input for the linear solver. The solution sol is initiated to 0, and then calculated by repeating the solve step. This process is called "projection".
k = -0.5
For tX = 1 To N
For tY = 1 To N
rhs(tX, tY) = k * (u(tX + 1, tY) - u(tX - 1, tY) + (v(tX, tY + 1) - v(tX, tY - 1)))
Next tY
Next tX
For tX = -1 To N2
For tY = -1 To N2
sol(tX, tY) = 0
Next tY
Next tX
For tX1 = 1 To 20
For tX = 1 To N
For tY = 1 To N
sol(tX, tY) = (rhs(tX, tY) + sol(tX - 1, tY) + sol(tX + 1, tY) +
sol(tX, tY - 1) + sol(tX, tY + 1)) / 4
Next tY
Next tX
Next tX1
For tX = 1 To N
For tY = 1 To N
u(tX, tY) = u(tX, tY) + k * (sol(tX + 1, tY) - sol(tX - 1, tY))
v(tX, tY) = v(tX, tY) + k * (sol(tX, tY + 1) - sol(tX, tY - 1))
Next tY
Next tX
Code Step 3
Yet further, a linear solver can be implemented by using a multigrid. If N1 is a power of 2, there are several levels of arrays rhs and sol in a multigrid, and each array is half the width and height of the previous. The coarsest level is solved over very few points and is interpolated onto the next, finer level. That one is then solved, and so on to the level with actual size N1xN1. It is an efficient algorithm. I implemented it in VB, thanks to Harald Köstler https://www10.cs.fau.de/publications/reports/TechRep_2008-03.pdf. VB is not the best language for a multigrid but this sample will work. The multigrid solver uses a calculation like this one on every level:
sol(tX, tY) = (rhs(tX, tY) + sol(tX - 1, tY) + sol(tX + 1, tY) +
sol(tX, tY - 1) + sol(tX, tY + 1)) / 4
It raises questions about the smoothness of the solution, but it creates curls on every level! Since these are interpolated to the next level, the result is a plane with larger and smaller curls.
To create curls in velocity, one or more of these three snippets can be used. They are similar, but not the same, so various effects can be generated.
Advection
We need yet another step for our fluid simulation. What we can see is the mixing of colours of fluid. The simple way is to add black (or grey) to a white fluid. The color array is usually called density. We start with a random distribution of density(dens) and of velocity (u and v). After playing with the velocity, as above, we calculate the moving of color or density according to that velocity. In fluid simulation, usually a process called advection is used. In a previous article, I already showed that backward advection can also produce curls. Here are some solutions on how to calculate the new distribution of density.
Code Step A, B
This is straightforward, no backward advection. Suited to investigate the effects above on their own.
For tX = 1 To N
For tY = 1 To N
tX1 = tX - u(tX, tY)
tY1 = tY - v(tX, tY)
dens(tX1, tY1) = dens0(tX, tY)
Next tY
Next tX
Code Step C, D
Just a little bit of backward advection. Above and here, u and v are floats, but tX, tY, tX1, tY1 are integers, the floats are rounded.
For tX = 1 To N
For tY = 1 To N
tX1 = tX - u(tX, tY)
tY1 = tY - v(tX, tY)
dens(tX, tY) = dens0(tX1, tY1)
Next tY
Next tX
Code Step E
Backward advection as it is used in most simulations. Because u and v are floats, the movement of density is usually not an integer shift. Hence, the interpolation over 4 points in between which the float distance falls. This produces smoother results than A to D above. The disadvantage is that some of the small curls get lost. A to D may be a bit rough around the edges and this one is more blurred.
N5 = N + 0.5
dt0 = dt
For tX = 1 To N
For tY = 1 To N
x = tX - dt0 * u0(tX, tY)
y = tY - dt0 * v0(tX, tY)
If x > N5 Then x = N5
If x < 0.5 Then x = 0.5
tX1 = Int(x)
tX2 = tX1 + 1
If y > N5 Then y = N5
If y < 0.5 Then y = 0.5
tY1 = Int(y)
tY2 = tY1 + 1
s1 = x - tX1
s0 = 1 - s1
t1 = y - tY1
t0 = 1 - t1
s2 = s0 * t0
s3 = s0 * t1
s4 = s1 * t0
s5 = s1 * t1
dens(tX, tY) = s2 * dens0(tX1, tY1) + s3 * dens0(tX1, tY2) +
s4 * dens0(tX2, tY1) + s5 * dens0(tX2, tY2)
u(tX, tY) = s2 * u0(tX1, tY1) + s3 * u0(tX1, tY2) + s4 * u0(tX2, tY1) + s5 * u0(tX2, tY2)
v(tX, tY) = s2 * v0(tX1, tY1) + s3 * v0(tX1, tY2) + s4 * v0(tX2, tY1) + s5 * v0(tX2, tY2)
Next tY
Next tX
This advection can be even more sophisticated. That requires a lot of processor time. Therefore, I chose not to include them here.
Code Step F
Forward advection with interpolation, I found it here.
For tX = 1 To N
For tY = 1 To N
sx1 = tX + u0(tX, tY)
sy1 = tY + v0(tX, tY)
If sx1 < 0 Then sx1 = 0
If sx1 > N1 Then sx1 = N1
If sy1 < 0 Then sy1 = 0
If sy1 > N1 Then sy1 = N1
rast(tX, tY, 1) = sx1
rast(tX, tY, 2) = sy1
Next tY
Next tX
For tX = 1 To N - 1
For tY = 1 To N - 1
rast(tX, tY, 3) = 0
Next tY
Next tX
For tX = 1 To N - 1
For tY = 1 To N - 1
sx1 = rast(tX, tY, 1)
sy1 = rast(tX, tY, 2)
sx2 = rast(tX + 1, tY, 1)
sy2 = rast(tX + 1, tY, 2)
sx3 = rast(tX + 1, tY + 1, 1)
sy3 = rast(tX + 1, tY + 1, 2)
sx4 = rast(tX, tY + 1, 1)
sy4 = rast(tX, tY + 1, 2)
mins = sx1
If sx2 < mins Then mins = sx2
If sx3 < mins Then mins = sx3
If sx4 < mins Then mins = sx4
tX1 = Int(mins)
mins = sy1
If sy2 < mins Then mins = sy2
If sy3 < mins Then mins = sy3
If sy4 < mins Then mins = sy4
tY1 = Int(mins)
mins = sx1
If sx2 > mins Then mins = sx2
If sx3 > mins Then mins = sx3
If sx4 > mins Then mins = sx4
tX2 = Int(mins)
mins = sy1
If sy2 > mins Then mins = sy2
If sy3 > mins Then mins = sy3
If sy4 > mins Then mins = sy4
tY2 = Int(mins)
a1 = sx1
a2 = sx2 - sx1
a3 = sx4 - sx1
a4 = sx1 - sx2 + sx3 - sx4
b1 = sy1
b2 = sy2 - sy1
b3 = sy4 - sy1
b4 = sy1 - sy2 + sy3 - sy4
aa = a4 * b3 - a3 * b4
For tX3 = tX1 To tX2
For tY3 = tY1 To tY2
bb = a4 * b1 - a1 * b4 + a2 * b3 - a3 * b2 + tX3 * b4 - tY3 * a4
cc = a2 * b1 - a1 * b2 + tX3 * b2 - tY3 * a2
det = bb * bb - 4 * aa * cc
If det >= 0 Then
det2 = Sqr(det)
If aa = 0 Then
w = 0
DoEvents
Else
w = (-bb + det2) / (2 * aa)
End If
wd = a2 + a4 * w
If wd <> 0 Then
z = (tX3 - a1 - a3 * w) / wd
If w >= 0 And w <= 1 And z >= 0 And z <= 1 _
Then
w1 = 1 - w
z1 = 1 - z
s2 = w1 * z1
s3 = w1 * z
s4 = w * z1
s5 = w * z
da = s2 * dens0(tX, tY) + s3 * dens0(tX + 1, tY) + _
s4 * dens0(tX, tY + 1) + s5 * dens0(tX + 1, tY + 1)
ua = s2 * u0(tX, tY) + s3 * u0(tX + 1, tY) + _
s4 * u0(tX, tY + 1) + s5 * u0(tX + 1, tY + 1)
va = s2 * v0(tX, tY) + s3 * v0(tX + 1, tY) + _
s4 * v0(tX, tY + 1) + s5 * v0(tX + 1, tY + 1)
dens(tX3, tY3) = (da + dens(tX3, tY3) * _
rast(tX3, tY3, 3)) / (rast(tX3, tY3, 3) + 1)
u(tX3, tY3) = (ua + u(tX3, tY3) * rast(tX3, tY3, 3)) / _
(rast(tX3, tY3, 3) + 1)
v(tX3, tY3) = (va + v(tX3, tY3) * rast(tX3, tY3, 3)) / _
(rast(tX3, tY3, 3) + 1)
rast(tX3, tY3, 3) = rast(tX3, tY3, 3) + 1
End If
End If
End If
Next tY3
Next tX3
DoEvents
Next tY
Next tX
That's it! For code V1, steps 1, 2 and 3 can be omitted, used on their own, or combined, of steps A, B, C, D, E, F, you take just one. V1 is random but smooth diffused density and velocity all over. V2 is specles of random density and velocity. Option G here is different from the above, especially to show the effect of steps 1, 2 and 3 on velocity. See the little hourglasses after 1 iteration? That shows the distribution of these steps, not so evenly diffused at all!
V3 is an initial sinusoid density that develops into the Rayleigh-Tailor instability. And it has nothing to do with solving a Navier-Stokes equation. (the formula for fluid flow):
I used step 1, that I invented myself to show the basic idea of the so-called "projection" step.
Playing with these steps, I have seen many different outcomes, some more like plants or crystals, some more like an abstract painting. What I was looking for, a texture like marbled stone or indeed stirred fluid, is, for example, the combination v1, 3, A.
Notes
- Run the demo from your file explorer, not from the IDE, to get "
App.Path" to point to the demo's directory. - Some patience when running the code! VB6 is slow but a good language for teaching purposes, I think.
- I got the suggestion to also calculate the movement of velocity in steps A or C. In step E and F, you can see that velocity moves as well. When I included it in A or C, the result is ugly IMHO, but I show them as B and D.
- In fluid simulation algorithms, there is often a diffuse step included. I left that out. Diffusion blurs a bit, and steps 1, 2, 3 and E, F do blur as well.
- To get pretty colored pictures, we assign a color to each of the possible density values (0 to 255).
- There is usually a lot of attention for boundaries. If you don't deal with them, the sides of your picture come caving in. What I do is define the arrays from -1 to N1 + 2, whereas the calculations are made between 1 and N1.
Conclusion
I think that real physics of fluid differs from these tricks to make velocity curl. It is a fact however that fluid velocity does curl, along with it moving straight forward. Either spinning or rolling: like a bullet in flight spins, or a maple seed that falls; or perhaps rolling like tumbleweed in the wind or wheels of a train on the track. Like when we say that waves roll in : fluid rolls. So to end here, a fluid-like warp (V4). Quite something else, but maybe more close to the physics of fluid surfaces than many fluid simulations suggest.
History
- 7th June, 2021: Initial version
