Solving Navier Stokes with Stable Fluids Algorithm in Python

#1: Discourse on Deep Learning Techniques for Engineering Physics

1. Introduction

Fluid mechanics is used as the standard mathematical framework on which these simulations are based. There is a consensus among scientists that the Navier-Stokes equations are a very good model for fluid flow.

In this article, I’ll discuss a reliable technique that fully resolves the Navier-Stokes equations for the first time, proposed by Jos Stam in his groundbreaking paper Stable Fluids [1]. This approach enables real-time interaction with three-dimensional fluids on a graphics workstation and is relatively simple to implement. This is accomplished by employing time steps that are significantly greater than those employed by Foster and Metaxas[2]. However, to obtain a stable solver Stable Fluid Algorithm (SFA) diverges from Foster and Metaxas’ approach of solution in order to build an unconditionally stable solver. SFA employs both Lagrangian and implicit techniques to solve the Navier-Stokes equations in place of Foster and Metaxas’ explicit Eulerian schemes. SFA was developed specifically for computer graphics applications; hence it cannot be found in the literature on computational fluids.

I’ll first discuss the theory behind solving Navier-Stokes using SFA, and then we will simulate hot fumes rising and swirling in an enclosed domain using the PhiFlow library in Python.

Once we have solved NS with SFA, we will use the methodology to optimize the solvers with Deep Learning techniques.

2. Navier-Stokes Equations

The Navier-Stokes equations are presented in this section, along with the adjustments that produced SFA’s stable solver. A velocity field (u) and pressure field (p) are used to represent a fluid whose density and temperature are virtually constant. These values often change throughout time and space depending on the boundaries that surround the fluid. We will use x to represent the spatial coordinate, which is x = (x,y) for two-dimensional fluids and x = (x,y,z) for three-dimensional fluids in both cases.

The Navier-Stokes equations [3] provide the evolution of these values over time, assuming that the velocity and pressure are known at some initial time t= 0:

The Navier-Stokes equations are derived by imposing that the fluid conserves both mass (Eq. 1) and momentum (Eq. 2). These equations must also be accompanied by boundary conditions. We will look at two relevant forms of boundary conditions: periodic boundary conditions and fixed boundary conditions. See Foster and Metaxas’ work for a good discussion of these boundary conditions in the case of a hot fluid [2]. In any case, the boundary conditions should be such that the normal component of the velocity field is zero at the boundary; no matter should traverse walls.

The Navier-Stokes equations show a relationship between the pressure and velocity fields. Combining Eq. 1 and 2 yields a single equation for velocity. We will look at how this single equation can be derived:

Any vector field w may uniquely be decomposed into the following form, according to a mathematical result known as the Helmholtz-Hodge Decomposition.

where u has zero divergences: ∇⋅u = 0;and q is a scalar field. Any vector field is the sum of a mass conserving field and a gradient field. This result allows us to define an operator P which projects any vector field w onto its divergence-free part u=Pw. The operator is, in fact, defined implicitly by multiplying both sides of Eq.3 by “∇”:

Neumann boundary condition is when the fluid lies in some bounded domain D. In such a case, the boundary conditions are given by a function defined on the boundary of the domain ∂D. Eq. 4 is a Poisson equation for the scalar field q with the Neumann boundary condition ∂q/∂n= 0 on ∂D. A solution to this equation is used to compute the projection u:

If we apply this projection operator on both sides of Eq. 2, we obtain a single equation for the velocity:

where we have used the fact that Pu = u and P∇p=0. This is the fundamental equation from which SFA is developed.

3. Stable Fluid Algorithm

Let us assume that the field has been resolved at a time t and that we wish to compute the field at a later time t+Δt. Eq. 5 is solved from an initial state of uₒ =u(x,0) by marching through time with a time step Δt in four steps:

Start:

From previous time step → wₒ(x) = u(x,t)

Step 1: Add Force

The easiest term to solve is the addition of the external force f (in Eq. 5). If we assume that the force does not vary considerably during the time step, then:

is a good approximation of the effect of the force on the field over the time step Δt. This is a good approximation in an interactive system since forces are only applied at the beginning of each time step.

Step 2: Advect

The term -(u ⋅ ∇)u in Eq. 5 describes how a disturbance in the fluid spreads from one place to another. This makes the Navier-Stokes equation non-linear.

In this step method of characteristics for solving partial differential equations is used. All fluid particles are transported by the fluid’s own velocity at each time step. As a result, we backtrace the point across the velocity field w₁over time Δt in order to determine the velocity at a point x at the new time t+Δt. This establishes a path p(x, s) that corresponds to a portion of the velocity field’s streamline. The particle’s current velocity at position x is then set to the velocity it had at its prior location at time Δt ago:

Step 3: Fourier Transform

Our algorithm takes a particularly simple form when we consider a domain with periodic boundary conditions. The periodicity allows us to transform the velocity into the Fourier domain:

In the Fourier domain, the gradient operator “∇” is equivalent to
the multiplication by ik, where i is the square root of -1.

Consequently, both the diffusion step (Step 4) and the projection step (Step 5) are much simpler to solve.

Step 4: Diffuse

The third step solves for the effect of viscosity and is equivalent to a diffusion equation:

To solve this, we discretize the diffusion operator ∇² and perform an Implicit time step.

Applying FFT:

Step 5: Project

The fourth step involves the projection step, which makes the resulting field divergence-free. From Eq. 4,

The projection step, therefore, requires a good Poisson solver. And we do this by applying Fourier Transforms.

Step 6: Inverse Fourier Transform

The solution at time t+Δt is then given by the last velocity field:

In fact, the SFA may be written down entirely in only a few lines. A quick Fourier transform and a particle tracer are all that are needed.

4. Moving Substances in a Fluid Domain

When a non-reactive material is introduced into the fluid, it will both diffuse and be advected by it. Examples of this phenomenon that are frequently seen include the patterns produced when coffee is mixed with cream or when cigarette smoke rises. Let a represent any scalar amount that is transported by the fluid. The temperature of a fluid, the density of dust, smoke, or cloud droplets, and a texture coordinate are a few examples of this quantity. An equation of the advection-diffusion type neatly describes the development of this scalar field:

where κₐ is a diffusion constant, αₐ is a dissipation rate, and Sₐ is
a source term. This equation is very similar in form to the Navier-Stokes equation. Indeed, it includes an advection term, a diffusion term, and a “force term” Sₐ. All these terms can be resolved exactly in the same way as the velocity of the fluid. The dissipation term not present in the Navier-Stokes equation is solved as follows over a time-step:

Now that we have discussed the theory and the algorithm. We will proceed to use PhiFlow to simulate a rising hot fume in an enclosed fluid domain.

5. Simulating Smoke Plume

The code below solves the incompressible Navier-Stokes equations in conjunction with an advection equation in a closed box.

Momentum: ∂u/∂t + (u ⋅ ∇) u = − 1/ρ ∇p + ν ∇²u + f

Incompressibility: ∇ ⋅ u = 0

Advection: ∂s/∂t + (u ⋅ ∇) s = α ∇²s + i

u: Velocity (2d vector)

p: Pressure

f: Forcing (here due to Buoyancy)

ν: Kinematic Viscosity

ρ: Density (here =1.0)

t: Time

∇: Nabla operator (defining nonlinear convection, gradient, and divergence)

∇²: Laplace Operators:

s: Concentration of a species (here hot smoke)

α: Diffusivity of the embedded species

i: Inflow of hot smoke into the domain.

Domain Setup:

• The domain is square and closed-off (wall BC everywhere)
• Initially, the fluid is at rest
• Initially, the concentration of smoke is zero everywhere
• There is a continuous inflow of hot smoke in a small circular patch in the bottom left of the domain
• The hot smoke exerts a force on the fluid due to Buoyancy
• This makes the fluid flow upwards and creates a plume pattern

Let's first install PhiFlow in a Jupyter Notebook setup.

Importing libraries:

Let’s define a few constants for the simulation.

We create grids for the quantities we want to simulate. For this example, we require a velocity field and a smoke density field.

The StaggeredGrid takes a couple of inputs. We first specify the initial velocity vectors at each grid cell which is 0 since the fluid is initially at rest. Next, we specify the extrapolation, which allows us to prescribe the boundary conditions, which are in turn modeled as the values of ghost cells around the domain. Since we are creating a fully closed domain, we specify the extrapolation to be ZERO which prescribes a wall boundary condition. We then set the extent of the domain in terms of discretization points 64×64 in the mesh grid. These grid points will be used on the physical domain, which is a 100×100 domain box.

Next, we define the grid for smoke which is our density species embedded into the fluid. Again the smoke will have 0 initial value. We will prescribe a Neumann bounday condition for the smoke by defining extrapolation as BOUNDARY. This center grid will be of high resolution with discretization points of 200×200. However, the domain extant will be the same as the staggered grid 100×100. This is because the grids for velocity and smoke will overlap.

The inflow will be used to inject smoke into a second centered grid smoke that represents the marker field from above.

The inflow sphere above is already using the “world” coordinates: it is located at x=25 along the first axis and y=15 on the second axis. The inflow field will have the same spatial dimension and resolution as the smoke field.

Let’s now create a step() function that will with our update routine. The step function takes the grid values for velocities and smoke from the previous step.

Step 1: Advect the density field (smoke) using the MacCormack method. Add inflow to prescribe “source” to the smoke field.

Step 2: Use the advected smoke field to compute buoyancy force via the Boussinesq model. The Boussinesq model uses multiplication by a tuple (0, buoyancy_factor) to turn the smoke field into a staggered, 2 component force field, sampled at the locations of the velocity components via the @ function. This function makes sure the individual force components are correctly interpolated for the velocity components of the staggered velocity. Note that this also directly ensures the boundary conditions of the original grid are kept. Put differently, we multiply smoke_next with (0.0, buoyancy_factor) which extends the smoke_next into a vector field with 0.0 in x component and buoyancy_force times the corresponding y scalar values in smoke_next. Then resample this vector field on the locations of the velocities in the velocity grid.

Step 3: We now advect the velocity according to the semi-Lagrangian scheme. Note: this is a self-advection term of NS equations. We then add the buoyancy force to the resulting velocities.

Step 4: We then add the diffusion term by implementing an explicit time step. At this point, the resulting velocities are of the next time step, but they are not incompressible yet.

Step 5: We now call the make_incompressible. This is typically the computationally most expensive step in the sequence above. It solves a Poisson equation for the boundary conditions of the domain and updates the velocity field with the gradient of the computed pressure.

The next cell takes one step and plots the density field after one simulation frame.

You’ll see that the hot smoke is appearing at the intended location. With this setup, we can easily advance the simulation forward in time a bit more by repeatedly calling the step function.

Let’s plot the density field at regular time intervals to see the evolution of the flow field.

Because of the inflow being located off-center to the left (with x position 30), the plume will curve towards the right when it hits the top wall of the domain. We can also take a look at the velocities. The steps listed above already store vector[0] and vector[1] components of the velocities as NumPy arrays, which we can show next.

6. Conclusion

PhiFlow is an incredibly intuitive framework for fluid simulations. In the near future, we will use PhiFlow to implement physics-embedded Deep Learning strategies. Till then, keep coding!

References

[1] Stam, Jos. (2001). Stable Fluids. ACM SIGGRAPH 99. 1999. 10.1145/311535.311548.

[2] N. Foster and D. Metaxas. Modeling the Motion of a Hot, Turbulent Gas. In Computer Graphics Proceedings, Annual Conference Series, 1997, pages 181–188, August 1997.

[3] A. J. Chorin and J. E. Marsden. A Mathematical Introduction to Fluid Mechanics. Springer-Verlag. Texts in Applied Mathematics 4. Second Edition., New York, 1990.