This package solves the inviscit, incompressible Navier Stokes equations in 2 or 3 (or any other) dimensions.

fluid_sim_demo.mp4

Installation

In a julia terminal, run the following commands to install this package:

using Pkg
Pkg.add(url="https://github.com/Alexander-Barth/FluidSimDemo.jl")

Then run the examples/fluid_sim_demo.jl.

Mathematical background

Navier-Stokes equations

The evolution equation of the velocity $\mathbf u$ is given by:

$$ \frac{d \mathbf u}{dt} + \mathbf u \cdot \nabla \mathbf u = -\frac{1}{\rho} \nabla p + \mathbf g $$

where $p$ is the pressure, $\rho$ is the constant density and $\mathbf g$ and external force like gravity. $\nabla$ is the nabla operator.

The flow is subjected to the following incompressibility constraint:

$$ \nabla \cdot \mathbf u = 0 $$

These equations are solved in this different steps:

1. Apply external forces $\mathbf g$ (in particular gravity)

$$ \mathbf {u'}^{(n)} = \mathbf u^{(n-1)} + \Delta t \mathbf g $$

Impose on boundaries $\mathbf {u'}^{(n)} = 0$ or the corresponding inflow/outflow velocity.

2. Advection of velocity field $\mathbf {u''}^{(n)}$. The combined force of inertia and external forces acting from $n-1$ to $n$ are written as $\mathbf F$:

$$ \mathbf {u''}^{(n)} = \mathbf {u}^{(n-1)} + \Delta t \mathbf F $$

3. The pressure is computed by requiring that the divergence of velocity must remain zero:

$$ \nabla \cdot \left( \mathbf F - \frac{1}{\rho} \nabla p^{(n)} \right) = 0 $$

which leads to:

$$ \nabla^2 p^{(n)} = \rho \nabla \cdot \mathbf {F}^{(n)} $$

Note that it is not necessary to compute $\mathbf F$ separatly as its divergence can be compute from $\mathbf {u''}^{(n)}$.

$$ \nabla \cdot \mathbf {F}^{(n)} = \frac{1}{\Delta t} \nabla \cdot \mathbf {u''}^{(n)} $$

$$ \nabla \cdot \mathbf {u''}^{(n)} = \Delta t \nabla \cdot \mathbf F $$

Choose pessure such that:

$$ \nabla \cdot \left( \mathbf {u''}^{(n)} - \Delta t \frac{1}{\rho} \nabla p^{(n)} \right) = 0 $$

$$ \nabla^2 p^{(n)} = \frac{\rho}{\Delta t} \nabla \cdot \mathbf {u''}^{(n)} $$

The pressure is solved iteratively (using Gauss-Seidel with overrelaxation) with a fixed number of iterations. In the following algorithm $\leftarrow$ is the assignement operator. For simplicity, the algorithm is outlined for the 2D case: Intitialize the iteration:

$$ \begin{alignat*}{2} & u'''_{i,j} &&\leftarrow u''^{(n)}_{i,j} \\ & v'''_{i,j} &&\leftarrow v''^{(n)}_{i,j} \\ & p_{i,j} &&\leftarrow 0 \\ \end{alignat*} $$

The time index $n$ is dropped as all parmeters are from the same time instance.

Compute pressure adjustement by:

$$ \Delta p \leftarrow -\frac{\rho \Delta x \Delta y}{4 \Delta t} \left(\frac{u'''_{i+1,j} - u'''_{i,j}}{\Delta x} + \frac{v'''_{i,j+1} - v'''_{i,j}}{\Delta y} \right) $$

Adjust the pressure

$$ p_{i,j} \leftarrow p_{i,j} + \Delta p $$

Update the velocity accordingly

$$\begin{array}{cc} u'''_{i,j} &\leftarrow u'''_{i,j} - \frac{\Delta p_{i,j} - \Delta p_{i-1,j}}{\Delta x} \frac{\Delta t}{\rho} \\ v'''_{i,j} &\leftarrow v'''_{i,j} - \frac{\Delta p_{i,j} - \Delta p_{i,j-1}}{\Delta y} \frac{\Delta t}{\rho} \\ \end{array} $$

After convergence, the velocity $u'''_{i,j}$ and $v'''_{i,j}$ correspond to the velocity for the next time step $\mathbf {u}^{(n)}$. On can show that $\mathbf {u}^{(n)}$ and $p$ satisfy the following equations:

$$ \mathbf {u}^{(n)} = \mathbf {u}^{(n-1)} + \Delta t \mathbf F - \frac{\Delta t}{\rho} \nabla p $$

Shallow-water model

The inviscid 2D shallow-water equations (derived from the Navier-Stokes equations) allow to compute the surface elevation η, and the u and v velocity components:

$$ \begin{align} \frac{\partial η}{\partial t} &+ \frac{\partial}{\partial x} \Bigl( (h + η) u \Bigr) + \frac{\partial}{\partial y} \Bigl( (h+η) v \Bigr) = 0,\\ \frac{\partial u}{\partial t} &+ u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} - f v = -g \frac{\partial η}{\partial x}\\ \frac{\partial v}{\partial t} &+ u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} + f u = -g \frac{\partial η}{\partial y} \end{align} $$

where h(x,y) is the depth of the water column, f is the Coriolis parameter and g is the acceleration due to gravity. Currently in the continuity equation (first equation), η (often in the order of 100 m) is neglected in front of η (order of cm).

Implementation

The julia code is implemented using a subset of julia that can be compiled to WebAssembly.

Credits

The 2D Navier Stokes is based on the compact implementation of Matthias Müller (which has been reimplemented in julia and extended in N-dimensions).