Multigrid applied to Incompressible Navier-Stokes Equation

Equations for fluid dynamics

In one of my research projects, I study numerical methods for fluid-structure interaction (FSI) problems. Examples of this type of problems are ubiquitous in nature, industries, and laboratories. For instance, as organisms swim in a liquid, the parts that generate propulsion, such as flagella and fins, are the solid structures that interact with the fluid environment around them. At the core of a numerical method for FSI problems, we need able to solve the PDEs that describe continuum mechanics. For fluid dynamics, we need to numerically solve the well-known Navier-Stokes equations. Assuming the fluid is incompressible Newtonian fluid, with constant density, we can write down two equations, first is the momentum balance equation, involving fluid velocity $\mathbf u$ , pressure $p$ , and body force acting on the fluid $\mathbf f$ :

$\frac{\partial \mathbf u}{\partial t} + (\mathbf u \cdot \nabla) \mathbf u= -\frac{1}{\rho}\nabla p + \nu \nabla^2 \mathbf u + \mathbf f$

where $\rho$ and $\nu$ are density and dynamic viscosity of the fluid. Second, we have the incompressibility constraint:

$\nabla \cdot \mathbf u = 0$

The first equation is a PDE in space and time, to solve it, we first discretize the spatial domain, then use finite difference to compute various spatial derivatives. Once the spatial derivatives are computed, we can use those to advance the solution forward in time. But the solution we get in this way is not guaranteed to fulfill the incompressibility constraint. To overcome this problem, we use Alexandre Chorin‘s projection method.

Projection method for incompressibility constraint

Schematic of leaving the admissible solution space of velocities, then project back on to it.

The mathematical formulation of the projection method would be a topic of another day, the idea is we would get an intermediate velocity update $\mathbf u^*$ by solving the incomplete momentum balance equation, now without the term involving pressure:

$\frac{\partial \mathbf u}{\partial t} + (\mathbf u \cdot \nabla) \mathbf u= \nu \nabla^2 \mathbf u + \mathbf f$

For example, if we were to use simple explicit time update rule, in discretized time step $n+1$ , we would have

$\frac{\mathbf u^* - \mathbf u^{n}} {\Delta t} = - (\mathbf u^n \cdot ) \mathbf u^n + \nu \nabla^2 {\mathbf u}^n + \mathbf f$

We update the velocities to $\mathbf u^*$ , then, we subtract in the pressure gradient term to get to next time step, i.e.

$\mathbf u^{n+1} = \mathbf u^* -\frac{\Delta t}{\rho} \nabla^2 p^n$

Now let’s apply divergence operator to both sides of the equation above, and have

$\nabla \cdot \left(\mathbf u^{n+1} - \mathbf u^* \right)=- \frac{\Delta t}{\rho} \nabla^2 p^n \implies$

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

The first term $\nabla \cdot \mathbf u^{n+1} =0$ because the incompressibility constraint needs to be enforced. Recognize this is a Poisson equation for pressure $p$ , with the source being

$\nabla \cdot \mathbf u^*$

We can choose how to would like to discretize the Laplacian and gradient operator, but in the end we would end up with a linear equation of the form

$\mathbf A \mathbf {\underline p} = \mathbf {\underline b}$

where $\mathbf {\underline p}$ and $\mathbf {\underline b}$ are discretized pressure and source term. This is a linear equation we can solve with multigrid.

In summary, we can use multigrid is in enforcing the incompressibility constraint in the Navier-Stokes equation via Chorin’s projection method. This method means we need to solve a Poisson equation at every time step. After discretization, the Poisson equation can be written in the form $\nabla\cdot \mathbf u^* = \frac{\Delta t}{\rho}\nabla^2 p^n$ . For details about multigrid algorithm, what it is and how it works, see the post Multigrid in a nutshell. In another post, I write about a project to parallelize the smoothing steps in multigrid using OpenMP and CUDA.

Yuexia Luna Lin

林月霞