Partial differential equations » Solving the Navier-Stokes equations

Steady incompressible Navier-Stokes with Picard iteration.

Introduction

The incompressible Navier-Stokes equations describe the motion of viscous fluids. In the steady-state case, they read:

where:

is the velocity field
is the pressure field
is the kinematic viscosity
is a body force
is the prescribed velocity on the boundary

The Navier-Stokes system is a mixed problem: we solve simultaneously for the velocity and pressure . The incompressibility constraint acts as a coupling condition.

The Stokes Problem

When the nonlinear convective term is dropped, we obtain the Stokes equations, which are linear:

The Stokes problem serves as the foundation for the Navier-Stokes solver and is also an important model in its own right (for very viscous or slow flows).

Weak Formulation

The weak formulation seeks such that for all :

where is the convective velocity (equal to at convergence). The mixed finite element spaces must satisfy the inf-sup (LBB) condition**. A classical choice is the Taylor-Hood pair:

Velocity: with polynomial order 2 (quadratic)
Pressure: with polynomial order 1 (linear)

The pressure is determined only up to a constant (only its gradient appears in the momentum equation). To fix the solution uniquely, we add a mean-pressure constraint** using a Lagrange multiplier from a P0 space.

Picard Iteration (Oseen Linearization)

The nonlinear term is handled by Picard iteration (also known as Oseen linearization). At each iteration :

Fix the convective velocity from the previous iteration
Solve the linearized (Oseen) problem for :
Repeat until convergence:

Implementation in Rodin

This example uses the Taylor-Green vortex as a manufactured solution. The code uses the H1 finite element space for higher-order elements.

Setup

#include <Rodin/Assembly.h>
#include <Rodin/Variational.h>
#include <Rodin/Variational/H1.h>
#include <Rodin/Solver/SparseLU.h>
#include <Rodin/IO/XDMF.h>

using namespace Rodin;
using namespace Rodin::Geometry;
using namespace Rodin::Variational;
using namespace Rodin::Solver;

Mesh and Spaces

We create a mesh on and define the Taylor-Hood finite element spaces. The H1 class supports arbitrary polynomial order via an std::integral_constant template parameter:

constexpr Real nu = 1.0;      // Kinematic viscosity
constexpr size_t maxIts = 12; // Max Picard iterations
constexpr Real tol = 1e-10;   // Convergence tolerance

const auto pi = Math::Constants::pi();

// Create mesh on unit square
Mesh mesh;
mesh = mesh.UniformGrid(Polytope::Type::Triangle, { 16, 16 });
mesh.scale(1.0 / 15.0);
mesh.getConnectivity().compute(1, 2);

// Taylor-Hood spaces: P2 for velocity, P1 for pressure
H1 uh(std::integral_constant<size_t, 2>{}, mesh, mesh.getSpaceDimension());
H1 ph(std::integral_constant<size_t, 1>{}, mesh);

// Lagrange multiplier for mean-pressure constraint
P0g p0g(mesh);

Manufactured Solution

The Taylor-Green vortex provides an exact solution for verification:

VectorFunction u_exact{
  sin(pi * F::x) * cos(pi * F::y),
  -cos(pi * F::x) * sin(pi * F::y)
};

auto p_exact = -0.25 * (cos(2 * pi * F::x) + cos(2 * pi * F::y));

VectorFunction f{
  2 * nu * pi * pi * sin(pi * F::x) * cos(pi * F::y)
    + pi * sin(2 * pi * F::x),
  -2 * nu * pi * pi * cos(pi * F::x) * sin(pi * F::y)
    + pi * sin(2 * pi * F::y)
};

Here F::x and F::y are symbolic coordinate functions that allow building expressions.

Trial and Test Functions

The Problem class supports multiple trial/test function pairs for mixed formulations. We define pairs for velocity, pressure, and the Lagrange multiplier:

TrialFunction u(uh);       // Velocity (vector)
TrialFunction p(ph);       // Pressure (scalar)
TrialFunction lambda(p0g); // Mean-pressure Lagrange multiplier

TestFunction  v(uh);       // Velocity test function
TestFunction  q(ph);       // Pressure test function
TestFunction  mu(p0g);     // Lagrange multiplier test function

Picard Iteration Loop

The convective velocity is stored as a GridFunction and updated after each solve. The convective term is computed using Mult (matrix-vector product) and Dot:

GridFunction w(uh);
w = VectorFunction{ Zero(), Zero() };  // Initial guess: zero velocity

GridFunction u_prev(uh);
u_prev = VectorFunction{ Zero(), Zero() };

for (size_t k = 0; k < maxIts; ++k)
{
  // Build the linearized (Oseen) problem
  Problem ns(u, p, lambda, v, q, mu);
  const auto conv_u = Mult(Jacobian(u), w);  // (w · ∇)u

  ns = nu * Integral(Jacobian(u), Jacobian(v))  // Viscous term
     + Integral(Dot(conv_u, v))                  // Convective term
     - Integral(p, Div(v))                       // Pressure-velocity coupling
     + Integral(Div(u), q)                       // Incompressibility
     + Integral(lambda, q)                       // Mean-pressure constraint
     + Integral(p, mu)                           // (symmetric counterpart)
     - Integral(f, v)                            // Body force
     + DirichletBC(u, u_exact);                  // Boundary conditions

  SparseLU solver(ns);
  solver.solve();

  // Check convergence
  P1 sh(mesh);
  GridFunction du(sh);
  du = Pow(Frobenius(u.getSolution() - u_prev), 2);
  Real updateNorm = std::sqrt(Integral(du).compute());

  // Update convective velocity for next iteration
  w.setData(u.getSolution().getData());
  u_prev.setData(u.getSolution().getData());

  if (updateNorm < tol)
    break;
}

Key Rodin expressions used:

Jacobian(u) — the gradient matrix
Mult(Jacobian(u), w) — the product , representing the convective term
Div(v) — the divergence
Frobenius(...) — the Frobenius norm for measuring convergence

XDMF Output

After convergence, export the velocity and pressure fields using the IO::XDMF class:

IO::XDMF xdmf("NavierStokes");
xdmf.grid().setMesh(mesh)
  .add("velocity", u.getSolution())
  .add("pressure", p.getSolution());
xdmf.write();

Open NavierStokes.xdmf in ParaView to visualize the velocity vectors and pressure contours.