Click here to Skip to main content
15,867,686 members
Articles / Programming Languages / Visual Basic 6
Tip/Trick

Fluid Rocks, Fluid Rolls

Rate me:
Please Sign up or sign in to vote.
5.00/5 (7 votes)
8 Jun 2021CPOL6 min read 4.1K   122   4   4
More on fluid simulation
In this tip, you will learn about the elements of many current fluid simulation algorithms that make the fluid velocity curl.

Combination V1, 1, F

Combination V1, 3, A

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:

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

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

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

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

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

VB.NET
N5 = N + 0.5
dt0 = dt ''* N
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.

VB.NET
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) ''quadrilateral with 4 corners
    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 ''integer bounds of quadrilateral
    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 ''multiply by matrix  [1 0 0 0; -1 1 0 0; -1 0 0 1; 1 -1 1 -1]
    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
        ''quadratic equation coeffs, aa*mm^2+bb*m+cc=0
      
            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) ''mapping of y on 0-1
                End If
                
                wd = a2 + a4 * w
                If wd <> 0 Then
                  z = (tX3 - a1 - a3 * w) / wd ''mapping of x on 0-1
                  If w >= 0 And w <= 1 And z >= 0 And z <= 1 _
                            Then ''tx3, ty3 is inside quadrilateral
                  
                    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!

Combination V2, 3, G

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):

Combination V3

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.

Combination V4

History

  • 7th June, 2021: Initial version

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)


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

Comments and Discussions

 
Questionsettings for picture 1 ? Pin
boobyy9-Sep-21 1:08
boobyy9-Sep-21 1:08 
AnswerRe: settings for picture 1 ? Pin
FranciskaR9-Sep-21 1:30
FranciskaR9-Sep-21 1:30 
PraiseRe: settings for picture 1 ? Pin
boobyy9-Sep-21 1:54
boobyy9-Sep-21 1:54 
GeneralMy vote of 5 Pin
LightTempler8-Jun-21 10:59
LightTempler8-Jun-21 10:59 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.