Approximation de problèmes mixtes

Model Problems

We consider now model problems as systems of PDEs where several functions are unknowns and which don’t play the same roles mathematically and physically.


\[ \left\{\begin{array}[c]{rl} -\Delta u + \nabla p & = f\ \mbox{ in } \Omega\\ \nabla \cdot u & = 0\ \mbox{ in } \Omega \end{array}\right. \] where \(u: \Omega \mapsto \RR^d\) is a velocity and \(p: \Omega \mapsto \RR\) is a pressure.


\[ \left\{\begin{array}[c]{rl} \sigma + \nabla u & = f\ \mbox{ in } \Omega\\ \nabla \cdot \sigma & = g\ \mbox{ in } \Omega \end{array}\right. \] where \(\sigma: \Omega \mapsto \RR^d\) is a velocity and \(u: \Omega \mapsto \RR\) is a hydraulic charge(pressure).


We shall focus on Stokes, but the abstract setting of the next section is the same for Stokes and Darcy.

Stokes and incompressible Navier-Stokes for Newtonian fluids

The Stokes model is the basis for fluid mechanics models and is a simplication of the Navier-Stokes equations where the viscous effects/terms are much bigger than the convective ones

\[ \left\{\begin{array}[c]{rl} \rho( \frac{\partial u}{\partial t} + u \cdot \nabla u) - \nu \Delta u + \nabla p & = f\ \mbox{ in } \Omega\\ \nabla \cdot u & = 0\ \mbox{ in } \Omega \end{array}\right. \] The first equation results from the conservation of momentum and the second from the conservation of mass.

The well-posedness of these problems results from a so-called which is not automatically transfered at the discrete level.

In practice In order to ensure that the finite element approximation is well-posed, we will need to choose approximation spaces that satisfy a compatibility condition that ensures that a discrete inf-sup condition is satisfied.

Saddle point problems

Abstract Continuous Setting


  • \(X\) and \(M\) two Hilbert spaces.[1]

  • two linear forms \(f \in X'=\mathcal{L}(X, \RR)\) and \(g \in M'=\mathcal{L}(M, \RR)\)

  • \(a \in \mathcal{L}(X\times X, \RR)\) and \(b \in \mathcal{L}(X\times M, \RR)\) two bilinear forms

We are interested in the following abstract problem:

Look for \((u,p) \in X \times M\) such that

\[ \left\{ \begin{array}[c]{rl} a(u,v) + b(v,p) & = f(v), \quad \forall v \in X\\ b(u,q) & = g(q), \quad \forall q \in M \end{array} \right. \]

Definition of a saddle point problem

If the bilinear form \(a\) is symmetric and positive on \(X\times X\), we say that [prob:chmixte:1] is a saddle point problem.

The structure of the problem is as follows

  • the space of solution is the same of the test space

  • the unknown \(p\) does not appear in the second equation

  • the unknown functions \(u\) and \(p\) are coupled via the same bilinear form \(b\) is the first and second equation.

The next question is :

Well posed problem


Let’s rewrite Problem [prob:chmixte:1].

Denote \(V=X\times M\) and introduce \(c \in \mathcal{L}(V\times V, \RR)\) such that

\[ cu,p),(v,q = a(u,v)+b(v,p)+b(u,q) \] and \(h\in \mathcal{L}(V,\RR)\) such that

\[ h(v,q) = f(v)+g(q) \] then problem [prob:chmixte:1] reads

Look for \((u,p) \in V\) such that

\[ \begin{array}[c]{rl} cu,p), (v,q & = h(v,q), \quad \forall (v,q) \in V \end{array} \]

We suppose that \(a\) is coercive over \(X\), the [prob:chmixte:2] is well-posed if and only if the bilinear form \(b\) satisfies the following inf-sup condition:

there exists \(\beta > 0\) such that

\[ \inf_{q \in M} \sup_{v \in X} \frac{b(v,q)}{||v||_X ||q||_M} \geq \beta \]

Lax-Milgram provides only a sufficient condition for well-posedness
The form \(c\) in [prob:chmixte:2] does not satisfy Lax-Milgram.

Let’s introduce the so-called Lagrangian \(l \in \mathcal{L}(X\times M, \RR)\) defined by

\[ l(v,q) = \frac{1}{2} a(v,v) + b(v,q) - f(v) - g(q) \]

We say that the point \((u,p)\in X\times M\) is a saddle point of \(l\) if

\[ \forall (v,q) \in X\times M, \quad l(u,q) \leq l(u,p) \leq l(v,p) \]

Under the hypothesys oF [thr:chmixte:1], the Lagrangian \(l\) defined by has a unique saddle point. Moreover this saddle point is the unique solution of problem [prob:chmixte:1].

Finite element approximation

Abstract Discrete Problem

We now turn to the approximation of the problem [prob:chmixte:1] by a standard Galerkin method in a conforming way.

Denote the two spaces \(X_h \subset X\) and \(M_h \subset M\), we consider the following problem:

Look for \((u_h,p_h) \in X_h \times M_h\) such that

\[ \left\{ \begin{array}[c]{rl} a(u_h,v_h) + b(v_h,p_h) & = f(v_h), \quad \forall v_h \in X_h\\ b(u_h,q_h) & = g(q_h), \quad \forall q_h \in M_h \end{array} \right. \]

We suppose that \(a\) is coercive over \(X\) and that \(X_h \subset X\) and \(M_h \subset M\).

Then the [prob:chmixte:3] is well-posed if and only if the following discrete inf-sup condition is satisfied:

there exists \(\beta_h > 0\) such that

\[ \inf_{q_h \in M_h} \sup_{v_h \in X_h} \frac{b(v_h,q_h)}{||v_h||{X_h} ||q_h||{M_h}} \geq \beta_h \]

The compatibility condition problem [prob:chmixte:3], to be well posed, requires that the spaces \(X_h\) and \(M_h\) satisfy the condition.

This is known as the Babuska-Brezzi (BB) or Ladyhenskaya-Babuska-Brezzi (LBB).

Regarding error analysis, we have the following lemma

Thanks to the Lemma of Céa applied to Saddle-Point Problems, the unique solution \((u,p)\) of problem [prob:chmixte:3] satisfies

\[ \begin{array}[c]{rl} ||u-u_h||X & \leq c{1h} \inf_{v_h \in X_h} ||u-v_h||X + c{2} \inf_{q_h \in M_h} ||q-q_h||M\\ ||p-p_h||_X & \leq c{3h} \inf_{v_h \in X_h} ||u-v_h||X + c{4h} \inf_{q_h \in M_h} ||q-q_h||_M \end{array} \] where

  • \(c_{1h} = (1+\frac{||a||_{X,X}}{\alpha})(1+\frac{||b||_{X,M}}{\beta_h})\) with \(\alpha\) the coercivity constant of \(a\) over X.

  • \(c_{2} = \frac{||b||_{X,M}}{\alpha}\)

  • \(c_{3h} = c_{1h} \frac{||a||_{X,X}}{\beta_h}\), \(c_{4h} = 1+ \frac{||b||_{X,M}}{\beta_h}+\frac{||a||_{X,X}}{\beta_h}\)

The constants \(c_{1h}, c_{3h}, c_{4h}\) are as large as \(\beta_h\) is small.

Linear system associated

The discretisation process leads to a linear system.

We denote

  • \(N_u = \dim {X_h}\)

  • \(N_p = \dim {M_h}\)

  • \(\{\phi_i\}_{i=1,...,N_u}\) a basis of \(X_h\)

  • \(\{\psi_k\}_{k=1,...,N_p}\) a basis of \(M_h\)

  • for all \(u_h = \sum_{i=1}^{N_u} u_i \phi_i\), we associate \(U \in \R{N_u}\), \(U=(u_1,\ldots,u_{N_u})^T\), the component vector of \(u_h\) is \(\{\phi_i\}_{i=1,\ldots,N_u}\)

  • for all \(p_h = \sum_{k=1}^{N_p} u_k \psi_k\), we associate \(P \in \R{N_p}\), \(P=(p_1,\ldots,p_{N_p})^T\), the component vector of \(p_h\) is \(\{\psi_k\}_{k=1,\ldots,N_p}\)

The matricial form of problem [prob:chmixte:3] reads

\[ \begin{bmatrix} \mathcal{A} & \mathcal{B}^T\\ \mathcal{B} & 0 \end{bmatrix} \begin{bmatrix} U \\ P \end{bmatrix} = \begin{bmatrix} F\\ G \end{bmatrix} \]

where the matrix \(\mathcal{A} \in \R{N_u,N_u}\) and \(\mathcal{B} \in \R{N_p,N_u}\) have the coefficients

\[ \mathcal{A}_{ij} = a(\phi_j,\phi_i), \quad \mathcal{B}_{ki} = b(\phi_i,\psi_k) \]

and the vectors \(\mathcal{F} \in \R{N_u}\) and \(\mathcal{G} \in \R{N_p}\) have the coefficients

  • \(F_i=f(\phi_i)\)

  • \(G_k=g(\psi_k)\)

  1. Since \(a\) is symmetric and coercive, \(\mathcal{A}\) is symmetric positive definite

  2. The matrix of the system is symmetric but not positive

  3. The inf-sup condition  is equivalent to the fact that \(\mathcal{B}\) is of maximum rank, i.e. \(\ker(\mathcal{B}^T) = \{0 \}\).

  4. From theorem [thr:chmixte:2], the matrix of the system  is invertible

When the inf-sup is not satisfied

The counter examples when the inf-sup condition  is not satisfied(e.g. \(\mathcal{B}\) is not maximum rank ) occur usually in two cases:


\(\dim {M_h} > \dim {X_h}\): the space of pressure is too large for the matrix \(\mathcal{B}\) to be maximum rank. In that case \(\mathcal{B}\) is injective (\(\ker(\mathcal{B}) = \{0\})\). We call this *locking*.

Spurious modes

there exists a vector \(Q^* \neq 0\) in \(\ker(\mathcal{B}^T)\). The discrete field\(q^*_h\) in \(M_h\), \(q^*_h=\sum_{k=1}^{N_p} Q^*_k \psi_k\), associated is called a *spurious mode*. \(q^*_H\) is such that

\[ b(v_h,q^*_h)=0. \]

We now introduce the Uzawa matrix as follows

The matrix

\[ \mathcal{U} = \mathcal{B} \mathcal{A}^{-1} \mathcal{B}^T \] is called the Uzawa matrix. It is symmetric positive definite from the properties of \(\mathcal{A}\), \(\mathcal{B}\)


The Uzawa matrix occurs when eliminating the velocity in system  and get a linear system on \(P\):

\[ \mathcal{U} P = \mathcal{B} \mathcal{A}^{-1} F - G \] then one application is to solve by solving iteratively and compute the velocity afterwards.

Mixed finite element for Stokes

Variational formulation

We start with the Well-posedness at the continuous level

  • We consider the model problem  with homogeneous Dirichlet condition on velocity \(u = 0\) on \(\partial \Omega\)

  • We suppose the \(f \in [L^2(\Omega)\)^d] and \(g \in L^2(\Omega)\) with \[ \int_\Omega g = 0 \] Introduce

    \[ L^2_0(\Omega) = \Big\{ q \in L^2(\Omega): \int_\Omega q = 0 \Big\} \]

The condition comes from the divergence theorem applied to the divergence equation and the fact that \(u=0\) on the boundary

\[ \int_\Omega g = \int_\Omega \nabla \cdot u = \int_{\partial \Omega} u \cdot n = 0 \] This is a necessary condition for the existence of a solution \((u,p)\) for the Stokes equations with these boundary conditions.

We turn now to the variational formulation.

The Stokes problem reads

Look for \((u,p) \in [H^1_0(\Omega)\)^d \times L^2_0(\Omega)] such that

\[ \left\{ \begin{array}[c]{rl} \int_\Omega \nabla u : \nabla v -\int_\Omega p \nabla \cdot v & = \int_\Omega f \cdot v, \quad \forall v \in [H1_0(\Omega)]d\\ \int_\Omega q \nabla \cdot u & = - \int_\Omega g q, \quad \forall q \in L^2_0(\Omega) \end{array} \right. \]

We retrieve the problem [prob:chmixte:1] with \(X=[H^1_0(\Omega)\)^d] and \(M=L^2_0(\Omega)\) and

\[ \begin{array}[c]{rlrl} a(u,v) &= \int_\Omega \nabla u : \nabla v,& \quad b(v,p) &= -\int_\Omega p \nabla \cdot v,\\ \quad f(v) &= \int_\Omega f \cdot v,& \quad g(q) &= - \int_\Omega g q \end{array} \]

Pressure up to a constant
The pressure is known up to a constant, that’s why we look for them in \(L^2_0(\Omega)\) to ensure uniqueness.

Finite element approximation

Denote \(X_h \subset [H^1_0(\Omega)\)^d] and \(M_h \subset L^2_0(\Omega)\)

Look for \((u_h,p_h) \in X_h \times M_h\) such that

\[ \left\{ \begin{array}[c]{rl} \int_\Omega \nabla u_h : \nabla v_h + \int_\Omega p_h \nabla \cdot v_h & = \int_\Omega f \cdot v_h, \quad \forall v_h \in X_h\\ \int_\Omega q_h \nabla \cdot u_h & = -\int_\Omega g q_h, \quad \forall q_h \in M_h \end{array} \right. \]

This problem, thanks to theorem [thr:chmixte:2] is well-posed if and only if \(X_h\) and \(M_h\) are such that there exists \(\beta_h > 0\)

\[ \inf_{q_h \in M_h} \sup_{v_h \in X_h} \frac{\int_\Omega q_h \nabla \cdot v_h}{||v_h||{X_h} ||q_h||{M_h}} \geq \beta_h \]

Some counter examples: bad finite element for Stokes

In this section, we present two classical bad finite element approximations.

Finite element \(\poly{P}_1/\poly{P}_0\): locking

Thanks to the Euler relations, we have

\[ \begin{array}[c]{rl} N_{\mathrm{cells}} - N_{\mathrm{edges}} + N_{vertices} &= 1-I\\ N^\partial_{\mathrm{vertices}} - N^\partial_{\mathrm{edges}} &= 0 \end{array} \]

where \(I\) is the number of holes in \(\Omega\).

We have that \(\dim {M_h} = N_{\mathrm{cells}}\),\(\dim {X_h} = 2 N^i_{\mathrm{vertices}}\) and so

\[ \dim {M_h} - \dim {X_h} = N_{\mathrm{cells}} - 2 N^i_{\mathrm{vertices}} = N^\partial_{\mathrm{edges}} - 2 > 0 \]

so \(M_h\) is too rich for the condition and we have \(\ker(\mathcal{B}) = \{0\}\) such that the only discrete \(u_h^*\), with components \(U^*\), satisfying \(\mathcal{B} U^*\) is the null field, \(U^*=0\).

Finite element \(\poly{Q}_1/\poly{P}_0\): spurious mode

We can construct in that case a function \(q_h^*\) on a uniform grid which is equal alternatively -1, +1 (chessboard) in the cells of the mesh, then

\[ \forall v_h \in [Q1_{c,h}]d, \quad \int_\Omega q^*_h \nabla \cdot v_h = 0 \] and thus the associated \(X_h\), \(M_h\) do not satisfy the condition.

Finite element \(\poly{P}_1/\poly{P}_1\): spurious mode

We can construct in that case a function \(q_h^*\) on a uniform grid which is equal alternatively -1, 0, +1 at the vertices of the mesh, then

\[ \forall v_h \in [P1_{c,h}]d, \quad \int_\Omega q^*_h \nabla \cdot v_h = 0 \] and thus the associated \(X_h\), \(M_h\) do not satisfy the condition.


The problem with the \(\poly{P}_1/\poly{P}_1\) mixed finite element is that the velocity is not rich enough.

A cure is to add a function \(v_h^*\) in the velocity approximation space to ensure that

\[ \int_\Omega q^_h \nabla \cdot v_h^ \neq 0 \] where \(q_h^*\) is the spurious mode.

To do that we add the bubble function to the \(\poly{P}_1\) velocity space.

Recall the construction of finite elements on a reference convex \(\hat{K}\). We say that \(\hat{b}: \hat{K} \mapsto \RR\) is a bubble function if:

  • \(\hat{b} \in H^1_0(\hat{K})\)

  • \(0 \leq \hat{b}(\hat{x}) \leq 1, \quad \forall \hat{x} \in \hat{K}\)

  • \(\hat{b}(\hat{C}) = 1, \quad \mbox{where} \hat{C}\) is the barycenter of \(\hat{K}\)


The function

\[ \hat{b} = (d+1)^{d+1} \Pi_{i=0}^d\ \hat{\lambda}_i \] where \((\hat{\lambda}_0, \ldots, \hat{\lambda}_d)\) denote the barycentric coordinates on \(\hat{K}\)

Denote now \(\hat{b}\) a bubble fonction on \(\hat{K}\), we set

\[ \hat{P} = [\poly{P}_1(\hat{K}) \oplus \mathrm{span} (\hat{b})]^d, \] and introduce

X_h &=& \Big\{ v_h \in [C^0(\bar{\Omega})]^d : \forall K \in \mathcal{T}_h, v_h
\circ T_K \in \hat{P}; v_{h_|{\partial \Omega}} = 0 \Big\}\\
M_h &=& P^1_{c,h}

The spaces \(X_h\) and \(M_h \cap L^2_0(\Omega)\) satisfy the compatibility condition  uniformly in \(h\).

Suppose that \((u,p)\), solution of [prob:chmixte:1], is smooth enough, ie. \(u \in [H^2(\Omega)\)^d \cap [H1_0(\Omega)]d] and \(p\in H^1(\Omega) \cap L^2_0(\Omega)\).

Then there exists a constant \(c\) such that for all \(h >0\)

\[ \| u- u_h \|{1,\Omega} + \|p-p_h\|{0,\Omega} \leq c h (\|u\|{2,\Omega} + \|p\|{1,\Omega}) \] and if the Stokes problem is stabilizing then

\[ \|u-u_h\|{0,\Omega} \leq c h^2 ( \|u\|{2,\Omega} +\|p\|_{1,\Omega}). \]

Stabilizing Stokes problem

We say that the Stokes problem is stabilizing if there exists a constant \(c_S\) such that for all \(f \in [L^2(\Omega)\)^d], the unique solution \((u,p)\) of with \(g=0\) is such that:

\[ \|u\|{2,\Omega} + \|p\|{1,\Omega} \leq c_S \|f\|_{0,\Omega} \] A sufficient condition for stabilizing Stokes problem is that the \(\Omega\) is a polygonal convex in 2D or of class \(C^1\) in \(\RR^d, d=2,3\).

Taylor-Hood Element

The mini-element solved the compatibility condition problem, but the error estimation in equation is not optimal in the sense that

  1. the pressure space is sufficiently rich to enable a \(h^2\) convergence in the pressure error,

  2. but the velocity space is not rich enough to ensure a \(h^2\) convergence in the velocity error.

The idea of the Taylor-Hood element is to enrich even more the velocity space to ensure optimal convergence in \(h\).

Here we will take \([\poly{P}_2\)^d] for the velocity and \(\poly{P}_1\) for the pressure.

Introduce \[\begin{aligned} \label{eq:chmixte:39} X_h &=& [P2_{c,h}]d\\ M_h &=& P^1_{c,h} \end{aligned} \]

The spaces \(X_h\) and \(M_h \cap L^2_0(\Omega)\) satisfy the compatibility condition  uniformly in \(h\).

Suppose that \((u,p)\), solution of problem [prob:chmixte:1], is smooth enough, ie. \(u \in [H^3(\Omega)\)^d \cap [H1_0(\Omega)]d] and \(p\in H^2(\Omega) \cap L^2_0(\Omega)\).

Then there exists a constant \(c\) such that for all \(h >0\)

\[ \| u- u_h \|{1,\Omega} + \|p-p_h\|{0,\Omega} \leq c h^2 (\|u\|{3,\Omega} + \|p\|{2,\Omega}) \] and if the Stokes problem is stabilizing then

\[ \|u-u_h\|{0,\Omega} \leq c h^3 ( \|u\|{3,\Omega} +\|p\|_{2,\Omega}). \]

Generalized Taylor-Hood element

We consider the mixed finite elements \(\poly{P}_k/\poly{P}_{k-1}\) and \(\poly{Q}_k/\poly{Q}_{k-1}\) which allows to approximate the velocity and pressure respectively with, on Simplices \[\begin{aligned} \label{eq:chmixte:42} X_h &=& [P{k}_{c,h}]d\\ M_h &=& P^{k-1}_{c,h} \end{aligned}\]] On Hypercubes \(\[\begin{aligned} \label{eq:chmixte:43} X_h &=& [Q^{k}_{c,h}\)^d\\ M_h &=& Q^{k-1}_{c,h} \end{aligned} \] We then have

\[ \|u-u_h\|{0,\Omega} + h ( \| u- u_h \|{1,\Omega} + \|p-p_h\|{0,\Omega} ) \leq c h^{k+1} (\|u\|{k+1,\Omega} +\|p\|_{k,\Omega}) \]

There are other stable discretization spaces

  • Discrete inf-sup condition: dictates the choice of spaces

  • Inf-sup stables spaces:

    • \(\mathbb Q_k\)-\(\mathbb Q_{k-2}\), \(\mathbb Q_k\)-\(\mathbb Q^{disc}_{k-2}\)

    • \(\mathbb P_k\)-\(\mathbb P_{k-1}\), \(\mathbb P_k\)-\(\mathbb P_{k-2}\), \(\mathbb P_k\)-\(\mathbb P^{disc}_{k-2}\)

    • Discrete inf-sup constant independent of \(h\), but dependent on \(k\)

Numerical validation: Test case

We consider the Kovasznay solution of the steady Stokes equations.

The exact solution reads as follows

\$\begin{array}{r c l} \mathbf{u}(x,y) & = & \left(1 - e^{\lambda x } \cos (2 \pi y), \frac{\lambda}{2 \pi} e^{\lambda x } \sin (2 \pi y)\right)^T \\ p(x,y) & = & -\frac{e^{2 \lambda x}}{2} \\ \lambda & = & \frac{1}{2 \nu} - \sqrt{\frac{1}{4\nu^2} + 4\pi^2}. \end{array}\$

The domain is defined as \$\domain = (-0.5,1) \times (-0.5,1.5)\$ and \$\nu = 0.035\$.

The forcing term for the momentum equation is obtained from the solution and is

\$ \mathbf{f} = \left( e^{\lambda x} \left( \left( \lambda^2 - 4\pi^2 \right) \nu \cos (2\pi y) - \lambda e^{\lambda x} \right), e^{\lambda x} \nu \sin (2 \pi y) (-\lambda^2 + 4 \pi^2) \right)^T\$

Dirichlet boundary conditions are derived from the exact solution.

1. An euclidian space which is complete for the norm induced by the scalar product