Ctrl+K

Stokes Equation

32. 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\)

32.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\).

32.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)
    gfu.vec.data = a.mat.Inverse(X.FreeDofs()) * f.vec
    return gfu

gfu = SolveStokes(X)
Draw (gfu.components[0], mesh, "velocity")
Draw (gfu.components[1], mesh, "pressure");

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, "velocity")
Draw (gfu.components[1], mesh, "pressure");   

---------------------------------------------------------------------------
NgException                               Traceback (most recent call last)
Cell In[5], line 4
      2 Q = L2(mesh, order=1)
      3 X = V*Q
----> 4 gfu = SolveStokes(X)
      5 Draw (gfu.components[0], mesh, "velocity")
      6 Draw (gfu.components[1], mesh, "pressure");   

Cell In[3], line 13, in SolveStokes(X)
     10 f = LinearForm( (x-0.5)*v[1]*dx).Assemble()
     12 gfu = GridFunction(X)
---> 13 gfu.vec.data = a.mat.Inverse(X.FreeDofs()) * f.vec
     14 return gfu

NgException: UmfpackInverse: Numeric factorization failed.

\(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, "velocity")
Draw (gfu.components[1], mesh, "pressure");

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, "velocity")
Draw (gfu.components[1], mesh, "pressure");

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).

Contents