Approximation de problèmes mixtes

Ce document est en anglais et doit être traduit.

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

1.1. 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 inf-sup condition 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.

2. Saddle point problems

2.1. 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:

Problem Mixed 1

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

2.1.1. Definition of a saddle point problem

Definition

If the bilinear form \(a\) is symmetric and positive on \(X\times X\), we say that the mixed problem 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 :

2.2. Well posedness

2.2.1. Reformulation

Let’s rewrite the mixed problem 1.

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

\[c((u,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 Problem Mixed 1 reads

Problem

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

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

Theorem

We suppose that \(a\) is coercive over \(X\), the Problem 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 Problem 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)\]

Definition

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

Proposition

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 Problem Mixed 1.

3. Finite element approximation

3.1. Abstract Discrete Problem

We now turn to the approximation of the Problem Mixed 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:

Formulation of the Abstract Discrete 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.\]

Theorem

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

Then the Formulation of the Abstract Discrete Problem 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 Formulation of the Abstract Discrete Problem, 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

lemma

Thanks to the Lemma of Céa applied to Saddle-Point Problems, the unique solution \((u,p)\) of problem Formulation of the Abstract Discrete Problem 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} ||p-q_h||_M\\ ||p-p_h||_M & \leq c_{3h} \inf_{v_h \in X_h} ||u-v_h||_X + c_{4h} \inf_{q_h \in M_h} ||p-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.

3.2. 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 Formulation of the Abstract Discrete Problem 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)\)

Since \(a\) is symmetric and coercive, \(\mathcal{A}\) is symmetric positive definite
The matrix of the system is symmetric but not positive
The inf-sup condition is equivalent to the fact that \(\mathcal{B}\) is of maximum rank, i.e. \(\ker(\mathcal{B}^T) = \{0 \}\).
From theorem Theorem, 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

Définition: Matrice d’Uzawa

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.

4. Mixed finite element for Stokes

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

Problem

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 [H^1_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 recover the formulation of Problem Mixed 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.

4.2. Finite element approximation

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

Problem

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 Theorem 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\]

4.3. Bad finite elements for Stokes

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

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

4.3.2. 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 [Q^1_{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.

4.3.3. 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 [P^1_{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.

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

Definition: Mini-Element

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

\[\begin{aligned} 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} \end{aligned}\]

lemma

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

Theorem

Suppose that \((u,p)\), solution of Problem Mixed 1, is smooth enough, ie. \(u \in [H^2(\Omega)]^d \cap [H_0^1(\Omega)]^d\) and \(p\in H^1(\Omega) \cap L_0^2(\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}).\]

Definition: 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\).

4.5. Taylor-Hood Element

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

the pressure space is sufficiently rich to enable a \(h^2\) convergence in the pressure error,
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 &=& [P^2_{c,h}]^d\\ M_h &=& P^1_{c,h} \end{aligned}\]

lemma

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

Theorem

Suppose that \((u,p)\), solution of problem Problem Mixed 1, is smooth enough, ie. \(u \in [H^3(\Omega)]^d \cap [H_0^1(\Omega)]^d\) and \(p\in H^2(\Omega) \cap L_0^2(\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\)

5. Test Cases for Stokes

5.1. Kovasznay

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 manufactured from the exact solution.

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