Thread by @jose_morval

It's time to solve some fluid dynamics equations! These are the well-known Navier-Stokes Equations which describe the motion of a viscous fluid. #ue5 #niagara #gamedev #techart #houdini #math

https://twitter.com/jose_morval/status/1580283593015730176

There are three or four unknows in this system of equations, depending on whether we are in 2D or 3D case. The main unknows are: the velocity (two or three components) and the pressure.

The first equations, one for each component of the velocity, are called "momentum equations" and are basically Newton's F=ma applied to a little fluid element. The last equation is the "continuity equation" and it means that the flow is incompressible, i.e., density constant.

The "velocity" is a vector field that represents how fast and in which direction the fluid moves through the space. The "pressure" is related with the force exerted by the fluid at any point inside it.

These equations need one more thing to "close" the system: a boundary condition. Usually, we don't solve the fluid equations in all (infinite) space. We restrict the equations to a finite region.

Given that our equations may vary in time and space we need to say something about the equations at some initial time (at every point in the space) and what happen at the "boundary" (in every moment).

For viscous fluids this reduces to: all velocity components equal zero at boundary region and initial velocity arrows inside our region (maybe zero too).

When we solve numerically this system of equations we obtain simulations like the one in the gif (thanks @keenanisalive).

This is all the fluid dynamics knowledge we need to know before we start. There are infinite books/blogs/papers about this topic but I would hightlight Joe Stam's Stable Fluids paper d2f99xq7vri1nk.cloudfront.net/legacy_app_fil… and the NVIDIA GPU Gems blog developer.nvidia.com/gpugems/gpugem…

Ok, now (really) we are going to solve (numerically) the Navier-Stokes equations. But...how? What means "solve" in this context? Let's start with something simpler.

I'm going to introduce two "simplifications" of the Navier-Stokes: one is Laplace-Poisson equation and the other is the heat equation. It can be observed, that both equations can be obtained when some of the terms of the Navier-Stokes equations vanish.

Solve numerically these equations consists in finding a finite set of values that (almost) satisfy the equations. The most common strategy is discretize the region of values that varies (i.e. space; and time, in heat eqn.).

Math says that in "the limit" of subdivision of our region, the numerical solution "converges" to the "real" or continuum one.

In the picture we have a visualization of "discretization". All figures represent the same function, u(x,y) = x^2-y^2 (that satisfy Laplace eqn.), with differents "levels" os discretization.

But here, we are only taking a known function and evaluating in our discrete region. Let's think in the other way: what conditions must satisfy our discrete points to solve Laplace eqn? Well, basically the expression in the picture.

If you have some math background, you'd know that these derivatives can be aproximated using some finite differences, so the last relation can be written as the picture below (we are taking f=0). We use "h" as the uniform distance between (x.y) points in the grid.

This is an interesting expression. It tells us that the value of a function that satisfy Laplace eqn. in a specific point is the average of the function evaluated in the adjacent points. This is called "Mean Value Property" of Laplace eqn.

Doing the same thing for the heat equation (using finite differences for time derivative) we get the expression in the picture below. In this case we use "Delta x" to mean the same as "h" and "Delta t" the subdivision in time dimension.

So, we have numerical "schemes" for our problem! In the second one is easy to see that we can "leave alone" the value with uppper script "n+1", that is the value of the function in (i,j) point in space and (n+1) in time. We use the "old" (time n) values to reach the next.

But, in the scheme for Laplace equation, we don't have time. The values at the left and the right of the equal sign have the same priority. For these cases, we have "iterative methods" like Jacobi or Gauss-Seidel that try to reach the solution bit a bit.

The scheme using this method reads as follows.

We are going to prototype this iterative method in Houdini (suited for nice visualization). We'll use u(x,y)=x^2-y^2 in the boundary to compare our last example. The pictures show the Laplace eqn implementation.

Here the solver over different size grids.

The solver for heat equation works (very) similarly. All the things we are been doing in this thread is the "core" of the fluid dynamics simulation that we are going to implement in Niagara, as we'll see soon.

Sorry for the long thread mixing, maybe, too hard or too easy concepts. But I think is good to see this kind of things (maybe in messy way) and take some time to think about how and why the things work. I think is good point to stop.

PD: In the visualization of the solver we are seeing the solver (with different grid sizes) versus its analytical solution u(x,y)=x^2-y^2. That is, evaluate the (x,y) point in the functions directly.