Partial differential equations » Solving the Navier-Stokes equations

Transient nonlinear lid-driven cavity with PETSc SNES.

Image

Introduction

The incompressible Navier-Stokes equations describe viscous flow with a nonlinear convective term. The PETSc example solves a transient lid-driven cavity problem on the unit square:

\[ \left\{ \begin{aligned} \rho\partial_t\mathbf{u} + \rho(\mathbf{u}\cdot\nabla)\mathbf{u} -\mu\Delta\mathbf{u} + \nabla p &= 0 && \mathrm{in}\ \Omega, \\ \nabla\cdot\mathbf{u} &= 0 && \mathrm{in}\ \Omega, \\ \mathbf{u} &= \mathbf{g} && \mathrm{on}\ \partial\Omega. \end{aligned} \right. \]

The top boundary moves with a prescribed horizontal lid velocity and the remaining walls use no-slip boundary conditions. Velocity is discretized with quadratic H1 elements and pressure with linear H1 elements, giving the Taylor-Hood pair.

Weak Formulation

Backward Euler gives the nonlinear residual at every time step:

\[ R(\mathbf{u}, p; \mathbf{v}, q) = { \rho \over \Delta t } \int_\Omega (\mathbf{u}^{n+1} - \mathbf{u}^n)\cdot\mathbf{v}\,dx + \rho\int_\Omega ((\mathbf{u}^{n+1}\cdot\nabla)\mathbf{u}^{n+1})\cdot\mathbf{v}\,dx + \mu\int_\Omega \nabla\mathbf{u}:\nabla\mathbf{v}\,dx - \int_\Omega p\,\nabla\cdot\mathbf{v}\,dx + \int_\Omega q\,\nabla\cdot\mathbf{u}\,dx + \epsilon\int_\Omega p q\,dx \]

The small pressure term $ \epsilon(p,q) $ fixes the pressure nullspace for the standalone example. The nonlinear state is initialized with the previous time step before each SNES solve.

PETSc SNES

PETSc SNES drives the nonlinear solve. The example assembles the Newton Jacobian explicitly, including both convective derivative terms:

\[ \rho\int_\Omega \left(((\mathbf{u}\cdot\nabla)\delta\mathbf{u}) + ((\delta\mathbf{u}\cdot\nabla)\mathbf{u}))\right) \cdot\mathbf{v}\,dx. \]

In Rodin notation these terms are:

const auto stateConvection = Mult(Jacobian(uState), uState);
const auto newtonConvection =
  Mult(Jacobian(du), uState) + Mult(Jacobian(uState), du);

Problem navierStokes(du, dp, v, q);
navierStokes =
    (rho / dt) * Integral(du, v)
  + rho * Integral(Dot(newtonConvection, v))
  + mu * Integral(Jacobian(du), Jacobian(v))
  - Integral(dp, Div(v))
  + Integral(Div(du), q)
  + epsP * Integral(dp, q)
  + (rho / dt) * Integral(uState, v)
  - (rho / dt) * Integral(uOld, v)
  + rho * Integral(Dot(stateConvection, v))
  + mu * Integral(Jacobian(uState), Jacobian(v))
  - Integral(pState, Div(v))
  + Integral(Div(uState), q)
  + epsP * Integral(pState, q)
  + DirichletBC(du, -uState).on(walls)
  + DirichletBC(du, lidVelocity - uState).on(top);

The SNES wrapper synchronizes PETSc's current nonlinear iterate into uState and pState before each residual and Jacobian assembly.

Running

Build the PETSc examples and run:

./examples/PETSc/PDEs/PETSc_Seq_NavierStokes \
  -ns_n 32 -ns_nt 100 -ns_T 1.0 \
  -snes_monitor -snes_converged_reason \
  -ksp_type preonly -pc_type lu

The example writes:

  • NavierStokes.mesh
  • NavierStokes_velocity.gf
  • NavierStokes_pressure.gf
  • NavierStokes.xdmf

Related Oseen Examples

The older semi-implicit examples now use Oseen names to reflect that they freeze the transport velocity and solve one linearized problem per step:

  • examples/PETSc/PDEs/Seq_OseenFlow.cpp
  • examples/PETSc/PDEs/Seq_OseenLidDrivenCavity.cpp

See Also

« Solving the Helmholtz equation | Partial differential equations | Poisson equation with periodic boundary conditions »