2D Stokes problem#

Keywords: PDE:Stokes Method:FEM FE:CrouzeixRaviart FE:TaylorHood Condition:Dirichlet Solver:Direct

We consider the problem of a flow entrained in a cavity, say \(\Omega=]0,1[\times]0,1[\):

\[\begin{split}\DeclareMathOperator{\divop}{div} \begin{cases} -\nu\bigtriangleup \mathbf{u} + \nabla p=\mathbf{f} & \mathrm{in\;}\Omega \\ \divop \mathbf{u} = 0 & \mathrm{in\;}\Omega \\ \mathbf{u}=\mathbf{g} & \mathrm{on\;}\Sigma=\{0\}\times ]0,1[ \\ \mathbf{u}=\mathbf{0} & \mathrm{on\;}\Gamma=\partial \Omega\backslash \Sigma \end{cases}\end{split}\]

with \(\mathbf{f} \in L^2(\Omega)\) and \(\mathbf{g}=(0,v_0),\ v_0\in H^{\frac{1}{2}}_{00}(\Sigma)\). Note that \(\mathbf{g}|\mathbf{n}=0\).

This problem has the following variational formulation: find \((\mathbf{u},p)\in H^1(\Omega)^2\times L^2_0(\Omega)\) such that for all \((\mathbf{v},p)\in H^1_0(\Omega)^2\times L^2_0(\Omega)\):

\[\begin{split}\DeclareMathOperator{\divop}{div} \begin{cases} \displaystyle \nu\int_\Omega \nabla \mathbf{u}:\nabla \mathbf{v} -\int_\Omega p \divop \mathbf{v}-\int_\Omega \divop \mathbf{u} q=\int_\Omega \mathbf{f} \cdot \mathbf{v} \\ \mathbf{u}=\mathbf{g} \mathrm{\;on\;}\Sigma \mathrm{\;and\;} \mathbf{u}=\mathbf{0} \mathrm{\;on\;}\Gamma \end{cases}\end{split}\]

The following XLiFE++ implementation of this problem solves the 2D Stokes problem using different couples of approximation for the velocity \(v\) and the pressure \(p\) : (P2, P0), (P2, P1), (CR, P0); CR corresponds to a P1 Crouzeix-Raviart finite element approximation leading to a non-conforming approximation of \(H^1\). Note that the (P1, P0) approximation is not working because of dof locking.

Using Paraview with the Glyph filter that represents 2D vectors with arrows, we get the following pictures:

Stokes

Fig. 11 Solution of the 2D Stokes problem using different approximations#

Note that to approximate the space \(L^2_0(\Omega)\), an essential condition of the zero mean type of pressure p has been used. Such a condition is quite costly because it interacts with all the pressure dofs. If we remove it, the matrix becomes singular (the pressure being defined to within a constant) and direct solvers fail. But iterative solvers do converge to the right velocity, and pressure to the right constant, determined by the initial guess of the iterative solver !

Note also that, when using Crouzeix-Raviart element, in order to visualize the velocity, it is required to project the CR solution in a standard Lagrange finite element space (here \(P_1\)).