Feel++

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.

Stokes

\[ \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.

Darcy

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

Applications

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

Denote

  • \(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

Reformulation

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:

Locking

\(\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}\)

Applications

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.

Mini-Element

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}\)

Example

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