Incompressible Navier-Stokes equations

The incompressible Navier-Stokes equations describe conservation of mass and conservation of momentum, which can be written as a divergence-free constraint and an evolution equation:

$$\begin{array}{rl}\mathrm{\nabla }\cdot u& =0,\\ \frac{\partial u}{\partial t}+\mathrm{\nabla }\cdot (u{u}^{\mathsf{T}})& =-\mathrm{\nabla }p+\nu {\mathrm{\nabla }}^{2}u+f,\end{array}$$

where $\mathrm{\Omega }\subset {\mathbb{R}}^d$ is the domain, $d\in \{2,3\}$ is the spatial dimension, $u=(u^1,\dots ,u^d)$ is the velocity field, $p$ is the pressure, $\nu$ is the kinematic viscosity, and $f=(f^1,\dots ,f^d)$ is the body force per unit of volume. The velocity, pressure, and body force are functions of the spatial coordinate $x=(x^1,\dots ,x^d)$ and time $t$. We assume that $\mathrm{\Omega }$ is a rectangular domain.

The equations are stated here in dimensionless form with a reference length and velocity of unity, in which case the viscosity is the inverse of the Reynolds number: $\nu =1/Re$. In the code, the viscosity is passed directly, e.g. params = (; viscosity = 1e-3).

Integral form

The integral form of the Navier-Stokes equations is used as starting point to develop a spatial discretization:

$$\begin{array}{rl}\frac{1}{|\mathcal{O}|}{\int }_{\partial \mathcal{O}}u\cdot n\mathrm{d}\mathrm{\Gamma }& =0,\\ \frac{\mathrm{d}}{\mathrm{d}t}\frac{1}{|\mathcal{O}|}{\int }_{\mathcal{O}}u\mathrm{d}\mathrm{\Omega }& =\frac{1}{|\mathcal{O}|}{\int }_{\partial \mathcal{O}}(-u{u}^{\mathsf{T}}-pI+\nu \mathrm{\nabla }u)\cdot n\mathrm{d}\mathrm{\Gamma }+\frac{1}{|\mathcal{O}|}{\int }_{\mathcal{O}}f\mathrm{d}\mathrm{\Omega },\end{array}$$

where $\mathcal{O}\subset \mathrm{\Omega }$ is an arbitrary control volume with boundary $\partial \mathcal{O}$, normal $n$, surface element $\mathrm{d}\mathrm{\Gamma }$, and volume size $|\mathcal{O}|$. We have divided by the control volume sizes in the integral form, so that all terms have the same units as their differential counterparts.

Boundary conditions

The boundary conditions on a part of the boundary $\mathrm{\Gamma }\subset \partial \mathrm{\Omega }$ are one or more of the following:

• Dirichlet: $u=u_{\text{BC}}$ on $\mathrm{\Gamma }$ for some $u_{\text{BC}}$;

• Neumann: $\mathrm{\nabla }u\cdot n=0$ on $\mathrm{\Gamma }$;

• Periodic: $u(x)=u(x+\tau )$ and $p(x)=p(x+\tau )$ for $x\in \mathrm{\Gamma }$, where $\mathrm{\Gamma }+\tau \subset \partial \mathrm{\Omega }$ is another part of the boundary and $\tau$ is a translation vector;

• Stress free: $\sigma \cdot n=0$ on $\mathrm{\Gamma }$, where $\sigma =(-p\mathrm{I}+2\nu S)$.

See Problem setup for how to prescribe boundary conditions in the code.

Pressure equation

Taking the divergence of the momentum equations yields a Poisson equation for the pressure:

$$-{\mathrm{\nabla }}^{2}p=\mathrm{\nabla }\cdot (\mathrm{\nabla }\cdot (u{u}^{\mathsf{T}}))-\mathrm{\nabla }\cdot f$$

Note the absence of time derivatives in the pressure equation. While the velocity field evolves in time, the pressure only changes such that the velocity stays divergence free.

If there are no pressure boundary conditions, the pressure is only unique up to a constant. Since only the gradient of the pressure appears in the equations, this constant can be set to zero without affecting the velocity field.

Other quantities of interest

Kinetic energy

The local and total kinetic energy are defined by $k=\frac{1}{2}\parallel u{\parallel }_2^2$ and $K=\frac{1}{2}\parallel u{\parallel }_{L^2(\mathrm{\Omega })}^2={\int }_{\mathrm{\Omega }}k\mathrm{d}\mathrm{\Omega }$. In the absence of viscosity, boundaries, and body forces, the total kinetic energy is conserved. The discretization used in this package preserves this property (see Spatial and temporal discretization).

Vorticity

The vorticity is defined as $\omega =\mathrm{\nabla }\times u$.

In 2D, it is a scalar field given by

$$\omega =-\frac{\partial u^1}{\partial x^2}+\frac{\partial u^2}{\partial x^1}.$$

In 3D, it is a vector field given by

$$\omega =\left(\begin{array}{c}-\frac{\partial u^2}{\partial x^3}+\frac{\partial u^3}{\partial x^2}\\ -\frac{\partial u^3}{\partial x^1}+\frac{\partial u^1}{\partial x^3}\\ -\frac{\partial u^1}{\partial x^2}+\frac{\partial u^2}{\partial x^1}\end{array}\right).$$

Note that the 2D vorticity is equal to the $x^3$-component of the 3D vorticity.