Stationary Stokes equations

Let \(\Omega \subseteq \mathbb{R}^d\), \(d = 2\) or \(d=3\), be the domain of the fluid and \(\partial \Omega\) its boundary. The boundary satisfies \(\partial \Omega = \Gamma_D \cup \Gamma_N\), where \(\Gamma_D\) denotes the portion of the boundary where the Dirichlet boundary conditions are imposed, and \(\Gamma_N\) the portion where Neumann boundary conditions are imposed. We note \(u : \Omega \mapsto \mathbb{R}^d\) and \(p : \Omega \mapsto \mathbb{R}\) respectively the velocity and the pressure of the fluid. Let \(\mu \in \mathbb{R}^{+}\) be its dynamic viscosity. The stationary Stokes equations link these variables by the following system:

where \(f : \Omega \mapsto \mathbb{R}^d\) represents external forces and \(D(u) = \frac{1}{2}(\nabla u + \nabla u^{T})\) is the deformation tensor.

Variational formulation

We suppose \(f \in L^2(\Omega)^d\) and choose \(v \in H^1(\Omega)^d\), \(q \in L^2_0(\Omega) = \{ q \in L^2(\Omega) | \int_{\Omega} q = 0\}\) as the test functions. The variational formulation of \ref{p:1} is given by :

Find \((u,p) \in H^1(\Omega)^d \times L^2_0(\Omega)\) such that:

We denote:

where \(a\) is a bilinear form and \(b\) is a linear form. The variational formulation can be rewritten as:

Let \(X = H^1(\Omega)^d, M = L^2_0(\Omega)\). Find \((u,p) \in X \times M \) such that:

Using the Lax-Milgram theorem, one can prove that this variational formulation has a unique solution \((u,p)\). We approximate this solution by the Galerkin method using a \(\mathbb{P}_2/\mathbb{P}_1\) finite element discretization. We note \(X_h \subset X\) and \(M_h \subset M\) two Lagrangian finite element spaces of dimension \(N_u\) and \(N_p\).

We will detail two different methods to solve this variational problem. The differences between the two methods get obvious when implementing them. The first method requires the definition of a composite function space \(X_h \times M_h\) which results in more implicit integrals calculations, while in the second method the function spaces and integrals are defined one by one.

First resolution method

A first method is to solve directly the discrete variational problem:

Find \((u_h,p_h) \in X_h \times M_h\) such that:

Second resolution method

The second method consists in using the matrix formulation of the variational problem. Let \((\phi)_i\) be the basis of \(X_h\) and \((\varphi)_i\) the basis of \(M_h\). We have:

The discrete variational formulation is rewritten as :

Taking \(v_h = (\phi)_j\) and \(q_h = (\varphi)_j\), one find the matrix formulation:

where:

We will detail the implementation of these two resolution methods for stationary Stokes equations.

A C++ file contains the implementation of the resolution of Stokes equations. For its execution we must define two other files, a geo file for the definition of the geometry and a cfg file, Stokes_config.cfg. The latter contains the configuration of the problem: the name of the geometry file and the definition of some function expressions.

The execution in \(2\)D and \(3\)D is respectively realized by the command line:

We will first show the different steps of the C++ file and then its configuration in the cfg file.

Mesh

We start by loading the mesh:

Function spaces

First resolution method

For the first resolution method, when considering directly the discrete variational formulation, we define a composite function space: \(\mathbb{P}_2/\mathbb{P}_1/\mathbb{P}_0\).

The definition of this composite function space is given by:

We can thus define the different elements:

The \(\mathbb{P}_0\) Lagrangian elements \(\lambda_1\) and \(\lambda_2\) are used to impose the two conditions :

Second resolution method

While considering the matrix formulation one defines separately the three function spaces:

The matrix formulation is given in the previous paragraph. We add the blocks \(C\) and \(C^T\):

where \(C = (c(\psi_i, \varphi_j))_{i,j} = (\int_{\Omega} \psi_i \varphi_j)_{i,j}\), and \((\psi)_i\) a basis of the \(\mathbb{P}_0\) Lagrangian finite element space.

We define the different blocks of this matrix and of the two vectors:

Expressions We initialize the expressions defined in the cfg file:

Integrals and resolution

We can then define the integrals.

First resolution method

For this first method, we define a bilinear form a and a linear form l, adding to them the respective integrals. The position of the blocks is managed implicitly by the composite function space. In the names of the mathematical functions used, for example grad and gradt, the t at the end indicates that it is applied to a trial element.

The resolution is then done by the command:

Second resolution method

For the second method each block of the matrix A and the vector F is defined separately. When defining the integrals it will be necessary to indicate manually the position of each block.

Export of results and computation of errors

To export the results and thus to be able to visualize them with Paraview, we add the following lines:

We also add the computation of errors:

The cfg file

The cfg file defines the geometry and function expressions. An example is given by:

To verify our implementation, we consider different test cases. For each case, we perform a study of the convergence of errors. Using the spaces of approximations \(\mathbb{P}_2\)/\(\mathbb{P}_{1}\), one must find a convergence in \(h^3\) for the velocity error in \(L^2\) norm and a convergence in \(h^2\) for the sum between the pressure error in \(L^2\) norm and the velocity error in \(H^1\) norm.

We choose \(\mu = 1\) for all test cases. For \(2\)D simulations, we use mesh steps:

and in \(3\)D:

Homogeneous Dirichlet boundary conditions 2D

First, we consider Stokes equations with homogeneous Dirichlet conditions defined on the unit square. The exact solutions are given by:

and

The errors and convergence orders for the different mesh sizes are given in the following table:

The table shows that the computed orders of convergence are respectively equal to \(3\) for the velocity error in \(L^2\) norm and equal to \(2\) for the sum between the velocity error in \(H^1\) norm and the pressure error in \(L^2\) norm. We have thus found the theoretical orders of convergence.

Non-homogeneous Dirichlet boundary conditions 2D

Second we consider as domain a unit square and apply non-homogeneous Dirichlet boundary conditions. We define the exact solutions by:

and

The non-homogenous boundary conditions have the exact solution \(u\).

The errors for the different mesh sizes are shown in the following table:

We find once again the theoretical orders of convergence, so order \(3\) for the velocity and order \(2\) for the sum.

Mixed boundary conditions 2D

Third we consider mixed boundary conditions on the rectangle [\(0,l\)]\(\times\)[\(0,h\)], where \(l = 5\) its length and \(h = 1\) its height. Non-homogenous Dirichlet conditions are applied on the left side of the rectangle, denoted by \(\Gamma_{inlet}\), non-homogenous Neumann conditions on the right side, \(\Gamma_{outlet}\), and homogenous Dirichlet conditions on the two other sides, \(\Gamma_{wall}\).

The exact solutions are given by:

and

where \(D = \frac{v_{max}}{0.5h(h-0.5h)}\) and \(v_{max} = 1\) the maximal velocity. For the non-homogenous Dirichlet boundary conditions on \(\Gamma_{inlet}\) we use the exact solution of \(u\) and the Neumann conditions on \(\Gamma_{outlet}\) are defined by:

This example is detailed in the Feel++ documentation. By displaying the solution obtained for pressure and velocity, for a mesh size \(h=0.1\), we obtain the same results:

As expected, the errors are equal to \(0\):

Non-homogeneous Dirichlet boundary conditions 3D

Finally, for the \(3\)D test case, we use the unit cube considering non-homogeneous Dirichlet boundary conditions. The exact solutions are given by:

and

To define the boundary conditions, we use the exact solution of \(u\).

The values of the computed errors are given by the following table:

The table shows that the orders of convergence are stabilizing around \(3\) and \(2\).