23.1. Stokes Equation

The Stokes equation is the basic model for incompressible fluids. While several terms of the advanced Navier-Stokes equations are skipped, the Stokes equation allows to study the difficulty of incompressibility constraints.

The unknown fields are

\[\begin{eqnarray*} \text{velocity} & & u : \Omega \rightarrow {\mathbb R}^d \\ \text{pressure} & & p : \Omega \rightarrow {\mathbb R} \end{eqnarray*}\]

The Stokes equation is a system of equations

\[\begin{split} \begin{array}{ccccl} -\Delta u & + & \nabla p & = & f \\ \operatorname{div} u & & & = & 0 \end{array} \end{split}\]

with a given force density \(f : \Omega \rightarrow {\mathbb R}^2\). The first equation models balance of momentum, the second one is the incompressibility constraint. In the simplest case, we set boundary conditions \(u = 0\) on \(\partial \Omega\)

23.1.1. Variational Formulation

The weak form is: Find \(u \in V := [H_0^1(\Omega)]^d\) and \(p \in Q:= L_2^0(\Omega)\) such that

\[\begin{split} \begin{array}{ccccll} \int \nabla u \nabla v & + & \int \operatorname{div} v \, p & = & \int f v & \forall \, v \in V \\ \int \operatorname{div} u \, q & & & = & 0 & \forall \, q \in Q \end{array} \end{split}\]

The pressure is defined up to a constant, and the divergence equation is satisfied for constant test-functions, anyway. Thus, we need the space \(L_2^0\), which are \(L_2\)-functions with zero mean value. For the beauty of the equation, we have substituted \(p\) by \(-p\).

23.1.2. Finite Element Spaces

We discretize the Stokes equation by different finite element spaces. For the velocity we use continuous, \(H^1\)-conforming elements. For the pressure, which is a field in \(L_2\), discontinuous finite elements are natural. But, alternatively, one may also use continuous elements.

We try different combinations of finite element spaces for \(V\) and \(Q\).

from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh(unit_square.GenerateMesh(maxh=0.05))

Use continuous elements of order \(k=2\) for each velocity component, and piece-wise constants for pressure:

V = VectorH1(mesh, order=2, dirichlet=".*")
Q = L2(mesh, order=0)
X = V * Q

def SolveStokes(X):
    u,p = X.TrialFunction()
    v,q = X.TestFunction()
    a = BilinearForm(X)
    a += (InnerProduct(grad(u),grad(v))+div(u)*q+div(v)*p)*dx
    # a += -1e-10*p*q*dx
    a.Assemble()
    a.mat[X.components[0].ndof,X.components[0].ndof]+=-1  # fixing pressure in one point

    f = LinearForm( (x-0.5)*v[1]*dx).Assemble()

    gfu = GridFunction(X)
    try:
        gfu.vec.data = a.mat.Inverse(X.FreeDofs()) * f.vec
    except:
        print ("matrix singular")
    return gfu

gfu = SolveStokes(X)
Draw (gfu.components[0], mesh, vectors={"grid_size" : 20 })
Draw (gfu.components[1], mesh);

Try a \(P^2 \times P^{1,dc}\) pairing:

V = VectorH1(mesh, order=2, dirichlet=".*")
Q = L2(mesh, order=1)
X = V*Q
gfu = SolveStokes(X)
Draw (gfu.components[0], mesh, vectors={"grid_size" : 20 })
Draw (gfu.components[1], mesh);

matrix singular

\(P^{2+} \times P^{1,dc}\), i.e. adding cubic bubbles:

V = VectorH1(mesh, order=2, dirichlet=".*")
V.SetOrder(TRIG,3)
V.Update()
Q = L2(mesh, order=1)
X = V*Q

gfu = SolveStokes(X)
Draw (gfu.components[0], mesh, vectors={"grid_size" : 20 })
Draw (gfu.components[1], mesh);

The Taylor Hood element \(P^2 \times P^1\) with continuous pressure:

V = VectorH1(mesh, order=2, dirichlet=".*")
Q = H1(mesh, order=1)
X = V*Q

gfu = SolveStokes(X)
Draw (gfu.components[0], mesh, vectors={"grid_size" : 20 })
Draw (gfu.components[1], mesh);

We see that the pressure space must not be too rich, otherwise we obtain a singular system (or a bad solutions, in particular for the pressure).

23.1.3. Stokes as an optimization problem

We consider the constrained minimization problem:

\[ \min_{u \in [H_0^1]^2 \atop \operatorname{div} u = 0} \tfrac{1}{2} \| \nabla u \|^2_{L_2} - (f,u)_{L_2} \]

We use the method by Lagrange and define the Lagrange function

\[ L(u,p) = \tfrac{1}{2} \| \nabla u \|^2_{L_2} - (f,u)_{L_2} + (\operatorname{div} u, p)_{L_2} \]

Taking the variation w.r.t. the \(u\)-variable:

\[ \partial_v L(u,p) = (\nabla u, \nabla v)_{L_2} - (f,v)_{L_2} + (\operatorname{div} v, p) \]

and the variation with respect to \(p\):

\[ \partial_q L(u,p) = (\operatorname{div} u, q) \]

Setting all derivatives to zero leads to the variational formulation of the Stokes equations.

23.1.3.1. Solvability theory

The existence (and uniqueness) of the \(u\)-component is easy to prove: The sub-space \(\{ u \in [H_0^1]^d : \operatorname{div} u = 0 \}\) is a Hilbert-space again, and thus Lax-Milgram provides the existence of the unique minimizer.

The existence of the Lagrange parameter is non-trivial. It requires to show that the divergence operator is surjective

\[ \operatorname{div} : [H_0^1]^d \rightarrow L_2^0 , \]

i.e.

\[ \forall \, g \in L_2^0 \; \exists u \in [H_0^1]^d : \operatorname{div} u = g \qquad \text{with} \qquad \| u \|_{H^1} \leq c \, \| g \|_{L_2} \]

In the context of variational formulations surjectivity is usually expressed by the \(\inf \sup\) condition: There exists \(\beta > 0\) such that

\[ \inf_{p \in L_2^0} \sup_{u \in [H_0^1]^d} \frac{ \int \operatorname{div} u \, p \, dx}{\| u \|_{H_0^1} \| p \|_{L_2} } \geq \beta \]

A corresponding discrete \(\inf \sup\)-condition must be satisfied by the combination of finite element spaces:

\[ \inf_{p_h \in Q_h} \sup_{u_h \in V_h} \frac{ \int \operatorname{div} u_h \, p_h \, dx}{\| u_h \|_{H_0^1} \| p_h \|_{L_2} } \geq \beta \]

The linear system has the matrix structure

\[\begin{split} \left( \begin{array}{cc} A & B^T \\ B \end{array} \right) \left( \begin{array}{c} u \\ p \end{array} \right) = \left( \begin{array}{c} f \\ 0 \end{array} \right) \end{split}\]

The discrete \(\inf \sup\) condition guarantees that the \(B\) matrix has full rank.