Feel++

Approximation de problèmes coercifs

Dans ce chapitre, on s’intéresse à l’approximation de problèmes coercifs:

  • le Laplacien et des variantes sur les conditions aux limites et le l’équation elle-même;

  • l’élasticité linéaire qui permet de modéliser de petites déformations d’un milieu continu déformable.

Le Laplacien

On s’instéresse dans cette section à l’approximation élément fini conforme du problème suivant:

Formulation forte du Laplacien

On cherche \$u\$ telle que

\[ \begin{split} -\Delta u &= f \mbox{ dans } \Omega\\ u &= 0 \mbox{ sur } \partial \Omega \end{split} \]

Cadre Mathematique

On suppose que \$f \in L^2(\Omega)\$.

La formulation faible du problème Formulation forte du Laplacien est la suivante:

Formulation faible pour des conditions de Dirichlet homogènes

On cherche \$u \in H^1_0(\Omega)\$ telle que \[\label{eq:65} \int_\Omega \nabla u \cdot \nabla v = \int_\Omega f\, v,\quad \forall v \in H^1_0(\Omega) \]

Problème bien posé

Introduisons

  • \$V = H^1_0(\Omega)\$ doté de la norme \$\|\cdot\|_{1,\Omega}\$ telle que \$\|v\|_{1,\Omega} = (\|v\|^2_{0,\Omega} + \|\nabla v\|^2_{0,\Omega})^{1/2}\$

  • la forme bilinéaire \$a \in \mathcal{L}(V \times V, \R)\$ telle que \$a(u,v) = \int_\Omega \nabla u \cdot \nabla v \$

  • la forme linéaire \$\ell \in \mathcal{L}(V, \R)\$ telle que \$l(v) = \int_\Omega f \nabla v \$

Le problème Formulation faible pour des conditions de Dirichlet homogènes s’écrit sous forme abstraite

On cherche \$u \in V\$ telle que \[ a(u, v) = \ell(v), \quad \forall v \in V \]

L’espace \$V\$ est un espace de Hilbert et les formes \$a\$ et \$\ell\$ sont continues sur \$V\times V\$ et \$V\$ respectivement.

Il ne reste plus qu’à vérifier si le problème est bien posé (existence d’une solution unique).

Pour cela on utilise démontre la coercivité de la forme bilinéaire \$a\$ sur \$V \equiv H^1_0(\Omega)\$.

Ceci se fait grâce au lemme suivant:

Soit \$\Omega\$ un ouvert borné de \$\R\$. Il existe une constante \$c_\Omega\$ (dépendente de \$\Omega\$ telle que \[ \forall v \in H^1_0(\Omega),\quad \|v\|{0,\Omega} \le c\Omega \|\nabla v\|_{0,\Omega} \]

\$c_\Omega\$ est homogène à une longeur et peut être interprétée comme une mesure caractéristique de \$\Omega\$.

Grâce à l’inégalité de Poincaré, on a le résultat suivant

La forme bilinéaire \$a\$ du problème Formulation faible pour des conditions de Dirichlet homogènes est coercive

On note tout d’abord que par l’inégalité de Poincaré et la définition de \$\|\cdot\|_{1,\Omega}\$

\[ \|v\|^2_{1,\Omega} \le (1 + c^2_\Omega) \|\nabla v\|^2_{0,\Omega} \] On en déduit que

\[ \forall v \in H^1_0(\Omega),\quad a(v,v) = \|\nabla v\|^2_{0,\Omega} \ge \frac{1}{1+c^2_\Omega} \|v\|^2_{1,\Omega} \]

Le Lemme de Lax-Milgram [thr:12] permet alors de conclure sur l’existence d’une solution unique pour le problème Formulation faible pour des conditions de Dirichlet homogènes.

Approximation conforme

On utilise une approximation conforme par éléments finis de Lagrange.

On considère \$\Omega\$ un polygone ou polyhèdre régulier de \$\RR^2\$ ou \$\RR^3\$ respectivement et un maillage \$\calTh = \{K_e\}_{e=1...\Ne}\$ de \$\Omega\$.

On considère un élément fini de référence \$(\hat{K},\hat{P},\hat{\Sigma})\$ tel que \$\Pk \subset \hat{P}\$ et \$k+1 > \frac{d}{2}\$, voir le théorème [thr:16].

On note

\[ L^k_{c,h} = \{ v_h \in C^0(\bar{\Omega}); \forall K \in \mathcal{T}_h,\ v_h \circ T_K \in \hat{P}\} \] où \$T_K\$ est la tranformation géométrique de \$\hat{K}\$ dans \$K\$.

  • Si on utilise \$\hat{P}=\Pk{k}\$ on a \$L^k_{c,h} = P^k_{c,h}\$.

  • Si on utilise \$\hat{P}=\Qk{k}\$ on a \$L^k_{c,h} = Q^k_{c,h}\$.

Afin de construire un espace d’approximation conforme \$V_h \subset V =H^1_0(\Omega)\$ on prend \[ \label{eq:71} V_h = L^k_{c,h} \cap H^1_0(\Omega) \] c’est à dire que les fonctions de \$V_h\$ satisfont les conditions aux limites outre le fait d’être dans \$L^k_{c,h}\$.

Le problème discret s’écrit alors

Trouver \$u_h \in V_h\$ telle que \[ a(u_h,v_h) = \ell(v_h),\, \forall v_h \in V_h \]

qui est bien posé (existence et unicité de \$u_h\$) car \$a\$ est coercive sur \$V\$ et que \$V_h \subset V\$.

On a le résultat suivant:

On suppose que \$u\$ solution de Formulation faible pour des conditions de Dirichlet homogènes est dans \$H^{k+1}(\Omega) \cap H^1_0(\Omega)\$, alors il existe une constante \$c_1\$ telle que pour tout \$h\$

\[ \|u-u_v\|{1,\Omega} \le c_1 h^k |u|{k+1,\Omega} \] et il existe une constante \$c_2\$ telle que pour tout \$h\$

\[ \|u-u_v\|{0,\Omega} \le c_2 h^{k+1} |u|{k+1,\Omega} \]

La preuve de ([eq:73]) est obtenue grâce au lemme de Cea ([eq:cea]) et du théorème d’interpolation ([eq:62]).

La preuve de ([eq:74]) est obtenue grâce au lemme d’Aubin-Nitsche qui permet d’affirmer qu’il existe une constante \$c\$ telle que

\[ \|u-u_h\|{0,\Omega} \leq c h |u-u_h|{1,\Omega} \] et donc que ([eq:74]) se déduit de ([eq:74]).

Implémentation avec {feelpp}

Avec {feelpp}, on ne construit pas explicitement l’espace \$V_h\$ mais \$L^k_{c,h}\$.

Le traitement des conditions aux limites de Dirichlet du problème ([eq:64]) peut être effectué de diverses façons, nous en verrons une.

Commencons par le maillage, dans un premier temps nous définissons le type du maillage contenant soit des éléments de type simplexe (segment,triangle, tetrahèdre) ou de type hypercube (segment, quadrangle, hexahèdre).

  // un maillage de simplexe dans $\R$ telle que la transformation
  // géométrique $T_K,\ K \in \calTh$  soit affine
  typedef Mesh<Simplex<d,1> > mesh_type;

  // un maillage d'hypercube dans $\R$ telle que la transformation
  // géométrique $T_K,\ K \in \calTh$  soit affine en chacune des variables
  // typedef Mesh<Hypercube<d,1> > mesh\_type;

  // generate the mesh associated to the unit square $[0,1]^2$ using triangles
  auto mesh = unitSquare();
Le mot clé auto permet de faire de l’inférence de type, pour plus de détails consultez la page C++11 de Wikipedia.

Ensuite nous pouvons définir l’espace \$L^k_{c,h}\$,

  // Vh est une structure de donnée allouée dynamiquement
  auto Vh = Pch<1>( mesh );
  // u est un élément de Vh
  auto u = Vh->element();
  // u est un autre élément de Vh
  auto u = Vh->element();

À présent, nous définissons les formes bilinéaires \$a\$ et \$\ell\$ qui sont respectivement des formes bilinéaires et linéaires.

  auto a = form2( _test=Vh, _trial=Vh ); (1)
  a = integrate( _range=elements(mesh), _expr=gradt(u)*trans(grad(v)) ); (2)

  auto l = form1( _test=Vh ); (3)
  l = integrate( _range=elements(mesh), _expr=f*id(v) ); (4)
1 \$a \in \mathcal{L}(V_h \times V_h,\ \RR)\$
2 \$a = \sum_{e=1...\Ne} \int_{K_e} \nabla \varphi_j \cdot \nabla \varphi_i,\quad i,j=1...,\dim{V_h}\$
3 \$\ell \in \mathcal{L}(V_h,\ \RR)\$
4 \$\ell = \sum_{e=1...\Ne} \int_{K_e} f \varphi_i,\quad i=1...,\dim{V_h}\$

Afin de traiter les conditions aux limites de Dirichlet homogènes, on peut utiliser le mot-clé on qui permet de les imposer de manière forte.

  a += on(_range=boundaryfaces(mesh), _element=u, _rhs=l, _expr=constant(0.) );
Le mot-clé constant permet de transformer une type numérique ( entier, flottant) en une expression utilisable par le langage de {feelpp}. Notez également l’opération += qui permet de rajouter le traitement des conditions de Dirichlet tout en gardant les contributions précédentes. L’opération = aurait d’abord remis à \$0\$ les entrées de la matrice associée à \$a\$.

Enfin nous pouvons résoudre le problème [prob:6]

  a.solve( _rhs=l, _solution=u );

Le listing complet

Conditions aux limites

Conditions aux limites de Dirichlet non homogène

On suppose toujours \$f \in L^2(\Omega)\$ et on se donne une fonction \$g \in C^{0,1}(\partial \Omega)\$

[\$g\$ est Lipschitzienne sur \$\partial \Omega\$].

On s’intéresse au problème suivant:

On cherche \$u : \Omega \rightarrow \RR\$ telle que

\[ \begin{split} -\Delta u &= f \mbox{ dans } \Omega\\ u &= g \mbox{ sur } \partial \Omega \end{split} \]

L’hypothèse \$g \in C^{0,1}(\partial \Omega)\$ permet d’affirmer qu’il existe \$u_g \in H^1(\Omega)\$ telle que \$u_{g_{|\partial \Omega}} = g\$.

On se ramène au problème avec conditions de Dirichlet homogène en faisant le change d’inconnue \$u_0=u-u_g\$ et on s’intéresse au problème suivant:

On cherche \$u_0 \in H^1_0(\Omega)\$ telle que

\[ a(u_0,v) = \ell(v) - a(u_g,v),\quad \forall v \in H^1_0(\Omega) \]

Ce problème est bien posé d’après Lax-Milgram, voir section précédente.

On suppose que \$u\$ solution de [prob:8] est dans \$H^{k+1}(\Omega) \cap H^1_0(\Omega)\$, alors il existe une constante \$c_1\$ telle que pour tout \$h\$

\[ \|u-u_v\|{1,\Omega} \le c_1 h^k |u|{k+1,\Omega} \] et il existe une constante \$c_2\$ telle que pour tout \$h\$

\[ \|u-u_v\|{0,\Omega} \le c_2 h^{k+1} |u|{k+1,\Omega} \]

Avec {feelpp}, les conditions Dirichlet non-homogènes sont traitées par exemple avec le mot-clé on.

Conditions de Dirichlet non homogènes
  auto g = sin(2*pi*Px() ); (1)
  (2)
  ...
  a += on( _range=boundaryfaces(mesh), _expr=g ); (3)
1 définition de la fonction, p.ex \$g=sin(2 \pi x)\$
2 définition de \$a\$
3 ajout des conditions de Dirichlet non-homogènes
Il n’y a pas besoin de rajouter le terme \$a(u_g,v)\$ au second membre \$\ell(v)\$, il est pris en compte automatiquement par on.

Voici le listing complet de l’exemple du laplacien avec conditions de Dirichlet non-homogène

Condition aux limites de Neumann

Étant donnés un réel \$\mu\$ strictement positif, \$f \in L^2(\Omega)\$ et \$g \in L^2(\partial \Omega)\$, on s’intéresse au problème suivant:

On cherche \$u : \Omega \rightarrow \RR\$ telle que

\[ \begin{split} -\Delta u + \mu u &= f, \mbox{ dans } \Omega\\ \partial_\Next u &= g, \mbox{ sur } \partial\Omega \end{split} \]

où \$\partial_\Next u = \nabla u \cdot \Next = \sum_{i=1}^d n_i \partial_i u\$ dénote la dérivée normale de \$u\$ avec \$\Next=(n_1,...,n_d) \in \RR^d\$ la normale extérieure unitaire en un point du bord de \$\Omega\$.

La formulation faible s’écrit

On cherche \$u \in H^1(\Omega)\$ telle que

\[ a( u, v ) = \ell(v),\ \forall v \in H^1(\Omega) \] avec

\[ a( u, v ) = \int_\Omega \nabla u \cdot \nabla v + \mu u v \] et

\[ \ell( v ) = \int_\Omega f v + \int_{\partial\Omega} g v \]

On a

\[ a(v, v) = \int_\Omega \nabla v \cdot \nabla v + \mu v v \ge \min(1,\mu) \int_\Omega \nabla v \cdot \nabla v + v v = \min(1,\mu) \|v\|_{1,\Omega}, \quad \forall v \in H^1(\Omega) \] ce qui nous permet d’affirmer que \$a\$ est coercive sur \$H^1(\Omega)\$ et que le problème [prob:13] est bien posé grâce à Lax-Milgram.

On peut approcher le problème [prob:13] par des éléments finis de Lagrange.

On utilise la formulation développée dans la section Approximation conforme

On cherche \$u_h \in P^k_{c,h}\$ tel que \[ a_h(u_h,w_h)=\ell(w_h), \quad \forall w_h \in P^k_{c,h} \]

Par rapport aux conditions de Dirichet, les conditions de Neumann sont directement (naturellement) traitées par la formulation faible. Elles ne sont pas directement imposées dans l’espace et les fonctions tests peuvent prendre des valeurs non-nulles au bord. Ces conditions sont traitées de manière approchée et non pas exacte.

L’analyse du problème [prob:13_1] est identique par le théorème [thr:17] aux mêmes estimations que dans le cas Dirichlet homogène.

Cas \$\lambda=0\$

Le cas \$\lambda=0\$ présente quelques difficultés techniques.

On a

On cherche \$u : \Omega \rightarrow \RR\$ telle que

\[ \begin{split} -\Delta u &= f, \mbox{ dans } \Omega\\ \partial_\Next u &= g, \mbox{ sur } \partial\Omega \end{split} \]

On observe ici qu’une condition nécessaire d’existence de solution est que
\$\int_\Omega f + \int_{\partial \Omega} g = 0\$

de par le théorème de la divergence.

la deuxième observation est que la solution de [prob:13_2] est connue à une constante additive près: si \$u\$ est solution et \$c\in \RR\$ alors \$u+c\$ est également solution.

Il convient alors de chercher la solution dans l’espace fonctionel suivant

H^1_*(\Omega) = \{ v \in H^1(\Omega):\quad \int_\Omega v = 0 \}

La formulation faible du problème [prob:13_2] s’écrit alors

On cherche \$u\in H^1_*(\Omega)\$, telle que \[ a(u,v) = \ell(v) \quad \forall v\in H^1_*(\Omega) \] avec \$a(u,v)=\int_\Omega \nabla u \cdot \nabla v\$.

Le caractère bien posé du problème [prob:13_3] résulte de la coercivité de \$a\$ sur \$H^1_*(\Omega)\$ qui elle-même résulte du Lemme suivante

Lemme d’inégalité de Poincaré-Wirtinger

Soit \$\Omega\$ un ouvert borné de \$\RR^d\$, il existe une constante \$C_\Omega\$ telle que \[ \|v\| \leq C_\Omega \|\nabla v\|{0,\Omega},\quad \forall v \in H^1*(\Omega) \]

Le problème [prob:13_3] peut être approché par des éléments finis de Lagrange ce qui conduit au problème discret suivant

On cherche \$u_h \in P^k_{c,h}\$ telle que a_h(u_h,v_h)=\ell(v_h) \quad \forall v_h\in P^k_{c,h}

L’espace d’approximation \$P^k{c,h}\$ n’est pas conforme dans \$H^1_*(\Omega)\$

le problème [prob:13_4] est singulier, il a une infinité de solution. L’une d’entre elle peut être approchée par une méthode itérative de type gradient conjugué. Notons \$u^*_h\$ cette solution alors la solution à moyenne nulle est

\$u_h=u^*_h-\frac{1}{|\Omega|}\int_\Omega u^*_h\$

Il s’agit donc d’effectuer a posteriori du calcul un post-processing pour se ramener à la solution à moyenne nulle. Dans {feelpp} on utilisera la fonction mean

// ...
// u est la solution du problème discret
a.solve(_rhs=l,_solution=u);

// calcul de la valeur moyenne de u
double meanu = mean(_range=elements(mesh),_expr=idv(u))(0,0);

// calcul de la solution à moyenne nulle u = u - meanu
u.add(-meanu);

Conditions aux limites de Robin

Étant donnés un réel \$\mu\$ strictement positif, \$f \in L^2(\Omega)\$ et \$g \in L^2(\partial \Omega)\$, on s’intéresse au problème suivant:

On cherche \$u : \Omega \rightarrow \RR\$ telle que

\[ \begin{split} -\Delta u &= f, \mbox{ dans } \Omega\\ \mu u + \partial_\Next u &= g, \mbox{ sur } \partial\Omega. \end{split} \]

La formulation faible s’écrit

On cherche \$u \in H^1(\Omega)\$ telle que \[\label{eq:84} a( u, v ) = \ell(v),\ \forall v \in H^1(\Omega) \] avec

\[ a( u, v ) = \int_\Omega \nabla u \cdot \nabla v + \int_{\partial \Omega} \mu u v \] et

\[ \ell( v ) = \int_\Omega f v + \int_{\partial\Omega} g v \]

On a

\[ \begin{split} a(v, v) & = \int_\Omega \nabla v \cdot \nabla v + \int_{\partial\Omega} \mu v v \\ & \geq \min(1,\mu)\left( \int_\Omega \nabla v \cdot \nabla v
\int_{\partial\Omega} v v\right) \\ &\geq \min(1,\mu) \|v\|_{1,\Omega}, \quad \forall v \in H^1(\Omega) \end{split} \]

La forme bilinéaire \$a\$ est donc coercive et le problème [prob:15] est bien posé grâce à Lax-Milgram.

On considère le problème discret suivant

On cherche \$u_h \in P^k_{c,h}(\Omega)\$ telle que

\[ a_h(u_h,v_h) = \ell(v_h)\quad \forall v_h \in P^k_{c,h} \]

Le problème est bien posé (\$P^k_{c,h} \subset H^1(\Omega)\$).

Comme pour le problème avec conditions de Neumann, les fonctions tests peuvent prendre des valeurs non nulles au bord. Les conditions de Robin(ou Fourier) ne sont satisfaites que de manière approchée et non pas exacte.

La convergence de \$u_h\$ est donnée par le théorème [thr:17].

Considérons \$\Omega=[0,1\$^2] et les données \$\mu=0.01\$, \$f=1\$ et \$g=0\$.

Advection-diffusion-réaction avec diffusion dominante

On s’intéresse au problème suivant:

On cherche \$u : \Omega \rightarrow \RR\$ telle que

\[ \begin{split} -\nabla \cdot ( \mathbf{\alpha} \nabla u ) + \mathbf{\beta} \cdot \nabla u + \mu u &= f \mbox{ dans } \Omega\\ \mathcal{B} u &= 0 \mbox{ sur } \partial \Omega\\ \end{split} \]

où \$\mathcal{B}\$ est l’opérateur prenant en compte les conditions aux limites. La variation sur les conditions aux limites de la section [sec:cond-aux-limit] s’appliquent.

  • \$\mathcal{B}u=u\$ condition de Dirichlet

  • \$\mathcal{B}u=(\mathbf{\alpha}\nabla{u}) \cdot n\$ condition de Neumann

  • \$\mathcal{B}u=\gamma u + (\mathbf{\alpha} \nabla u )\cdot n\$ condition de Robin

Les différents termes de l’opérateur ont une interprétation physique

  • \$-\nabla \cdot ( \mathbf{\alpha} \nabla u )\$ est un terme de diffusion,

  • \$\mathbf{\beta} \cdot \nabla u\$ est un terme d’advection,

  • \$\mu u\$ est un terme de réaction (production ou destruction).

Ce type d’équation est très fréquente en ingéniérie, biologie ou encore finance.

Transfert de chaleur

\$u\$ is the temperature, \$\alpha=\kappa \mathcal{I}\$ où \$\kappa\$ est la conductivité thermique, \$\beta\$ est le champ d’écoulement, \$\mu = 0\$ et \$f\$ est chaleur apportée exterieurement par unité de volume.

Advection Diffusion

\$u\$ est la concentration d’un soluté transportée dans un champ d’écoulement \$\beta\$. La matrice \$\alpha\$ modélise la diffusivité du soluté soit de la diffusion moléculaire soit du mélange turbulent du fluide transporteur. La production ou destruction par réaction chimique is pris en compte par le terme linéaire \$\mu u\$. Le second membre \$f\$ modélise les sources ou puits.

On suppose que \$\mathbf{\alpha} \in [L^{\infty}(\Omega)]^{d\times d}\$, \$\mathbf{\beta} \in [L^{\infty}(\Omega)]^d\$, \$\nabla \cdot \mathbf{\beta} \in L^{\infty}(\Omega)\$ et \$\mu \in L^\infty(\Omega)\$.

On suppose que l’opérateur \$\mathcal{L}\$ tel que \$\mathcal{L} u = -\nabla \cdot ( \alpha \nabla u ) + \beta \cdot \nabla u + \mu u\$ est elliptique au sens suivant:

L’opérateur \$\mathcal{L}\$ est elliptique si il existe une constante \$\alpha_0\$ telle que presque pour tout \$x\in\Omega\$

\[ \forall \xi=(\xi_1, \ldots , \xi_n)\in\RR^n,\quad { \sum_{i,j=1}^n \alpha_{ij}(x) \, \xi_i \, \xi_j \ge \alpha_0 \, \| \xi \|^2 } \]

Le Laplacien est dans la catégorie des opérateurs elliptiques, il correspond à \$\mathbf{\beta} = 0\$, \$\mu = 0\$ et \$\mathbf{\alpha}=\mathcal{I}_d\$ avec \$\mathcal{I}_d\$ la matrice identité de \$\RR^{d\times d}\$.

Condition de Dirichlet homogène

Nous imposons \$u=0\$ sur \$\partial \Omega\$.

Nous multiplions \$\mathcal{L} u = f\$ par une fonction test (suffisamment régulière) s’annulant au bord et on intègre sur \$\Omega\$. En intégrant par parties (formule de Green) nous avons

\$\int_\Omega \nabla \cdot (\mathbf{\alpha} \nabla u ) v = \int_\Omega (\mathbf{\alpha} \nabla u ) \cdot \nabla v -\int_{\partial \Omega} ((\sigma \nabla u)\cdot n) v \nabla u) v + \mu u v\$

ce qui donne

\$\int_\Omega (\mathbf{\alpha} \nabla u ) \cdot \nabla v + (\mathbf{\beta} \cdot \nabla u) v + \mu u v = \int_\Omega f v\$

Afin que les intégrales sur \$\Omega\$ aient un sens, nous demandons par exemple que \$u\$ et \$v\$ aient la régularité suivante:

\$u\in H^1(\Omega)\quad\text{et}\quad v\in H^1(\Omega).\$

Comme \$u\in H^1(\Omega)\$, le théorème [toto] implique que \$u\$ ait une trace sur le bord. Comme \$u_{|\partial \Omega} = 0\$, on cherche en fait \$u\$ dans \$H^1_0(\Omega)\$. Les fonctions tests sont également prises dans \$H^1_0(\Omega)\$ ce qui donne la formulation faible suivante

On cherche \$u \in H^1_0(\Omega)\$ telle que

\[ a(u,v) = \ell(v) \quad \forall v \in H^1_0(\Omega) \] avec

\[ a(u,v)=\int_\Omega (\mathbf{\alpha} \nabla u ) \cdot \nabla v + (\mathbf{\beta} \cdot \nabla u) v + \mu u v \] et

\[ \ell(v) = \int_\Omega f v \]

Nous avons le résultat suivant .Une solution du problème faible est solution du problème fort

Si \$u\$ résout [eq:89], alors \$\mathcal{L}u=f\$ presque partout dans \$\Omega\$ et \$u=0\$ presque partout sur \$\partial \Omega\$.

Soit \$\phi\in \mathcal{D}(\Omega)\$ et \$u\$ une solution de [prob:12]. On a \[ \begin{split} \left\langle-\nabla \cdot (\alpha \nabla u), \phi \right\rangle_{\mathcal{D}',\mathcal{D}} &= \left\langle (\alpha \nabla u), \nabla \phi \right\rangle_{\mathcal{D}',\mathcal{D}} = \int_\Omega \alpha \nabla u \cdot \nabla \phi \\ & = \int_\Omega (f - \beta \cdot \nabla u - \mu u ) \phi, \end{split} \] ce qui donne \$\langle\mathcal{L}u,\phi\rangle_{\mathcal{D}',\mathcal{D}} = \int_\Omega f \phi\$. Or \$\mathcal{D}(\Omega)\$ est dense dans \$L^2(\Omega)\$, nous avons donc que \$\mathcal{L}u=f\$ dans \$L^2(\Omega)\$. C’est à dire que \$\mathcal{L}u=f\$ presque partout dans \$\Omega\$. Enfin \$u=0\$ presque partout sur \$\partial \Omega\$ par définition de \$H^1_0(\Omega)\$ d’après le théorème [toto] \$\blacksquare\$

Conditions de Dirichlet Non homogène

Nous imposons \$u=0\$ sur \$\partial \Omega\$ où \$g : \partial \Omega \rightarrow \RR\$ est une fonction donnée.

Nous supposons que \$g\$ soit suffisamment régulière telle qu’il existe un relèvement \$u_g\$ de \$g\$ dans \$H^1(\Omega)\$, c’est à dire qu’il existe une fonction \$u_g\$ de \$H^1(\Omega)\$ telle que \$u_g=g\$ sur \$\partial \Omega\$.

On cherche \$u \in H^1(\Omega)\$ telle que \[ u=u_g+\phi,\quad \phi \in H^1_0(\Omega) \] où \$\phi\$ est solution de

\[ a(\phi,v) = \ell(v) - a(u_g,v) \quad \forall v \in H^1_0(\Omega) \] avec

\[ a(u,v)=\int_\Omega (\mathbf{\alpha} \nabla u ) \cdot \nabla v + (\mathbf{\beta} \cdot \nabla u) v + \mu u v \] et

\[ \ell(v) = \int_\Omega f v \]

Soit \$g\in H^{\frac{1}{2}}(\partial \Omega)\$, si \$u\$ résout [prob:12_1], alors \$\mathcal{L}u=f\$ presque partout dans \$\Omega\$ et \$u=g\$ presque partout sur \$\partial \Omega\$.

La preuve est similaire à celle de [prob:12].

Lorsque l’opérateur \$\mathcal{L}\$ est le Laplacien, [prob:12_1] est appelé un problème de Poisson.

Conditions de Neumann

Soit une fonction \$g:\partial \Omega \rightarrow \RR\$, nous voulons imposer \$\mathcal{B}u=(\alpha \nabla u)\cdot n = g\$ sur \$\partial \Omega\$.

Dans le cas où \$\alpha = \mathcal{I}\$, la condition de Neumann correspond à spécifier la dérivée normale de \$u\$ car \$\nabla u \cdot n = \partial_n u = \frac{\partial u}{\partial n}\$.

Nous procédons de la même façon que précédemment et en utilisant la condition de Neumann dans l’intégrale de surface [eq:adr-green], on obtient la formulation faible suivante:

On cherche \$u \in H^1(\Omega)\$ telle que

\[ a(u,v) = \ell(v) \quad \forall v \in H^1(\Omega) \] avec

\[ a(u,v)=\int_\Omega (\mathbf{\alpha} \nabla u ) \cdot \nabla v + (\mathbf{\beta} \cdot \nabla u) v + \mu u v \] et

\[ \ell(v) = \int_\Omega f v + \int_{\partial \Omega} g v \]

Soit \$g\in L^2(\partial \Omega)\$, si \$u\$ résout [prob:12_2], alors \$\mathcal{L}u=f\$ presque partout dans \$\Omega\$ et \$\alpha \nabla u \cdot n = g\$ presque partout sur \$\partial \Omega\$.

En prenant les fonctions tests dans \$\mathcal{D}(\Omega)\$ , on a immediatement que \$\mathcal{L}u = f\$ presque partout dans \$\Omega\$. Nous avons donc que \$-\nabla \cdot (\alpha \nabla u ) \in L^2(\Omega)\$. Le corollaire [B59] implique que \$\alpha \nabla u \cdot n \in H^{\frac{1}{2}}(\partial \Omega)' = H^{-\frac{1}{2}}(\partial \Omega)\$ du fait que \[ \forall \phi \in H^{\frac{1}{2}}(\partial \Omega), \quad \left\langle\alpha \nabla u \cdot n, \phi\right\rangle_{H{-\frac{1}{2}},H{\frac{1}{2}}} = \int_{\Omega} -\nabla \cdot (\alpha \cdot \nabla u) u_\phi + \int_\Omega \alpha \nabla u \cdot \nabla u_\phi \] où \$u_\phi \in H^1(\Omega)\$ est un relèvement de \$\phi\$ dans \$H^1(\Omega)\$. On a alors que [prob12_2] donne \[ \forall \phi \in H^{\frac{1}{2}}(\partial \Omega), \quad \left\langle\alpha \nabla u \cdot n, \phi\right\rangle_{H{-\frac{1}{2}},H{\frac{1}{2}}} = \int_{\partial \Omega} g \phi \] montrant que \$\alpha \nabla u \cdot n = g\$ dans \$H^{-\frac{1}{2}}(\partial \Omega)\$ et donc, par conséquent, dans \$L^2(\Omega)\$ du fait que \$g\$ soit dans \$L^2(\Omega)\$.\$\blacksquare\$

Conditions mixtes Dirichlet-Neumann

documentation en cours

Conditions de Robin

documentation en cours

Soient deux fonctions \$g,\gamma : \partial \Omega \rightarrow \RR\$, nous voulons imposer \$\gamma u + (\alpha \nabla u)\cdot n = g\$ sur \$\partial \Omega\$. En utilisant la relation sur l’intégrale de surface [eq:adr-green], nous avons la formulation faible suivante

On cherche \$u \in H^1(\Omega)\$ telle que

\[ a(u,v) + \int_{\partial \Omega} \gamma u v = \ell(v) \quad \forall v \in H^1(\Omega) \] avec

\[ a(u,v)=\int_\Omega (\mathbf{\alpha} \nabla u ) \cdot \nabla v + (\mathbf{\beta} \cdot \nabla u) v + \mu u v \] et

\[ \ell(v) = \int_\Omega f v + \int_{\partial \Omega} g v \]

Soit \$g\in L^2(\partial \Omega)\$ et \$\gamma \in L^\infty(\partial \Omega\$. Si \$u\$ résout [prob:12_4], alors \$\mathcal{L}u=f\$ presque partout dans \$\Omega\$ et \$\gamma u + \alpha \nabla u \cdot n = g\$ presque partout sur \$\partial \Omega\$.

On procède de manière similaire aux preuves précédentes. \$\blacksquare\$

Nous récapitulons dans la table suivante les différentes formulations

Table 1. Formulation faible pour différentes conditions aux limites
Problème Espace Forme bilinéaire Forme linéaire

Dirichlet homogène

\$H^1_0(\Omega)\$

\$a(u,v)\$

\$\int_\Omega f v\$

Neumann

\$H^1(\Omega)\$

\$a(u,v)\$

\$\int_\Omega f v + \int_{\partial \Omega} g v\$

Dirichlet Neumann

\$H^1_{\partial \Omega_D}(\Omega)\$

\$a(u,v)\$

\$\int_\Omega f v + \int_{\partial \Omega_N} g v\$

Robin (Fourier)

\$H^1(\Omega)\$

\$a(u,v) + \int_{\partial \Omega} \gamma uv\$

\$\int_\Omega f v + \int_{\partial \Omega} g v\$

Résumé

Les problèmes considérés précédent se mettent tous sous la forme suivante, sauf le problème de Dirichlet non homogène

\$\left\{\begin{array}{l} \text{ Chercher } u \in V \text{ telle que}\\ a(u,v)=f(v), \quad \forall v \in V \end{array}\right.\$

où \$V\$ est un espace de Hilbert satisfaisant

\$H^1_0(\Omega) \subset V \subset H^1(\Omega)\$

Dans le cas de conditions de Dirichlet non-homogènes, nous avons \$u\in H^1(\Omega)\$ et \$u=u_g+\phi\$ où \$u_g\$ est un relèvement de la donnée au bord \$g\$ et \$\phi\$ résout le problème générique avec des conditions de Dirichlet homogènes.

Conditions aux limites essentielles et naturelles

Il est important d’observer le traitement différent entre les conditions de Dirichlet et Neumann ou Robin conditions.

Les conditions de Dirichlet sont imposées explicitement dans l’espace fonctionnel où la solution est recherchée, et les fonctions de test disparaissent (i.e. \$v=0\$) sur la partie correspondante de la frontière. Pour cette raison, les conditions de Dirichlet sont souvent appelées conditions aux limites essentielles .

Les conditions de Neumann et Robin ne sont pas imposées par le cadre fonctionnel mais par la formulation faible elle-même. Le fait que les fonctions test ont des degrés de liberté sur la partie correspondante de la frontière est suffisant pour faire respecter les conditions limites en question. Pour cette raison, ces conditions sont souvent appelées conditions aux limites naturelles.

Coercivité

Cette section est en cours d’écriture

Soient \$f\in L^2(\Omega), stem:[\mathbf{\alpha} \in [L^{\infty}(\Omega)]^{d\times d}\$, \$\mathbf{\beta} \in [L^{\infty}(\Omega)]^d\$, \$\nabla \cdot \mathbf{\beta} \in L^{\infty}(\Omega)\$ et \$\mu \in L^\infty(\Omega)\$.

On note \$p = \essinf_{x \in \Omega} (\mu -\frac{1}{2} \nabla \cdot \mathbf{\beta})\$ et \$c_\Omega\$ est la constante de l’inégalité de Poincaré [B23].

(i) Les problèmes avec conditions de Dirichlet homogènes [prob12] et non-homogènes [prob12_2] sont bien posés si

\[ \alpha_0 > - \min( 0, \frac{\gamma}{c_\Omega} ) \]

(ii) Le problème avec condition de Neumann [prob12_1] est bien posé si

\[ p > 0\quad\text{et}\quad \essinf_{x \in \Omega}(\beta \cdot n ) \geq 0 \]

(iii) Le problème avec condition de Dirichlet-Neumann [prob12_3] est bien posé si [eq:92] est vérifiée et

\[ \mathrm{meas}(\partial \Omega_D) > 0\quad\text{et}\quad \partial \Omega^- = \{ x\in \partial \Omega ; (\beta \cdot n)(x) < 0\} \subset \partial \Omega_D \]

(iv) On pose \$q = \essinf_{x \in \Omega} (\gamma +\frac{1}{2} \mathbf{\beta}\cdot n)\$. Le problème avec conditions de Robin [prob12_4] non-homogènes [prob12_2] est posé si

\[ p \geq 0,\quad q \geq 0,\quad\text{et}\quad pq \neq 0. \]

Nous prouvons (i) et (iv).

Preuve de (i) En utilisant l’ellipticité de \$\mathcal{L}\$ et l’identité suivante(Formule de Green) \[ \int_\Omega u(\beta \cdot \nabla u) = -\frac{1}{2} \int_\Omega (\nabla \cdot \beta) u^2 + \frac{1}{2} \int_{\partial \Omega} (\beta \cdot n) u^2 \] alors on a \[ \forall u \in H^1_0(\Omega),\quad a(u,u) \geq \alpha_0|u|^2_{1,\Omega} + p \|u\|^2_{0,\Omega}. \] En posant \$c=\min(0,\frac{p}{c_\Omega}\$ et en utilisant l’inégalité de Poincaré [B53] on a \[ \forall u \in H^1_0(\Omega),\quad a(u,u) \geq \left(\alpha_0+\frac{c}{c_\Omega}\right)|u|^2_{1,\Omega} \geq \sigma \|u\|^2_{1,\Omega}. \] avec \$\sigma=\frac{c_\Omega(c_\Omega\alpha_0+c)}{1+c^2_\Omega}\$.

Cela montre que \$a\$ est coercive sur \$H^1_0(\Omega)\$.

Le caractère bien posé résulte alors du Lemme de Lax-Milgram pour les conditions de Dicirichlet homogènes et de la proposition [TODO] pour les conditions de Dirichlet non-homogènes.\$\square\$

Preuve de (iv) Notons \$a(u,v)=\int_\Omega \alpha \nabla u\cdot \nabla u + (\beta \cdot \nabla u) v + \mu u v + \int_{\partial \Omega} \gamma u v\$. Nous avons l’inégalité suivante: \[ \forall u \in H^1(\Omega),\quad a(u,u) \geq \alpha_0|u|^2_{1,\Omega} + p \|u\|^2_{0,\Omega} + q\|u\|^2_{0,\partial \Omega}. \]

Si \$p>0\$ et \$q\geq 0\$, la forme bilinéaire est coercive sur \$H^1(\Omega)\$ avec comme constante stem:[\sigma = \min(\alpha_0,p).

Si \$p\geq 0\$ et \$q > 0\$, la forme bilinéaire est coercive sur \$H^1(\Omega)\$ grâce au Lemme [TODO].

Dans les deux cas, le caractère bien posé est obtenu par le Lemme de Lax-Milgram.\$\square\$ Ceci termine les preuves de (i) et (iv). \$\blacksquare\$

La coercivité est garantie si \$\alpha_0\$ est suffisamment grand c’est à dire que si la diffusion est dominante.

Approximation conforme dans \$H^1\$

Cette section est en cours d’écriture

L’approximation élément fini est similaire à celle du Laplacian, de plus les variantes sur les conditions aux limites s’appliquent également: condition de Dirichlet non homogène, de Neumann ou de Robin.

Soit \$\Omega\$ un polyèdre dans \$\RR^d\$, soit \$\{\mathcal{T}_h\}\$ une famille de maillages de \$\Omega\$, et soit \$\{\hat{K}, \hat{P}, \hat{E}\}\$ un élément fini Lagrange de référence du degré \$k \geq 1\$.

Soit \$P^k_{c,h}\$ l’espace d’approximation conformé \$H^1\$ défini par

\$P^k_{c,h}=\{v \in C^0(\bar{\Omega}); \quad \forall K \in \mathcal{T}_h v_{|K} \in \mathbb{P}^k (K) \}\$

Pour obtenir une approximation conforme dans \$V\$, nous rajoutons les conditions aux limites, i.e, nous avons

\$V_h = P^k_{c,h} \cap V\$

Dans le cas de conditions de Dirichlet homogène cela donne

\$V_h=\{ v_h \in P^k_{c,h} ;\ v_h = 0 \mbox{ sur } \partial \Omega\}\$

Dans le cas Neumann et Robin, nous avons \$Vh=P^k_{c,h}\$. Enfin dans le cas Dirichlet-Neumann, nous avons

\$V_h=\{v_h \in P^k_c; \ v_h = 0 \mbox{ sur } \partial \Omega_D\}\$

Le problème discret s’ecrit

Trouver \$u_h\$ dans \$V_h\$ telle que \[ a(u_h,v_h)=\ell(v_h), \quad \forall v_h \in V_h \]

Nous cherchons à présent à estimer l’erreur \$u-u_h\$ en norme \$H^1\$ et \$L^2\$.

Estimation \$H^1\$

Soit \$\Omega\$ un polyèdre dans \$\RR^d\$, soit \$\{\mathcal{T}_h\}\$ une famille de maillages de \$\Omega\$, et soit \$\{\hat{K}, \hat{P}, \hat{E}\}\$ un élément fini Lagrange de référence du degré \$k \geq 1\$. Nous avons \$\lim_{h\rightarrow 0} \|u-u_h\|_{1,\Omega} = 0\$ et pour \$\frac{d}{2} < s < k+1, il existe une constante stem:[c\$ telle que \[ \forall h, \quad \|u-u_h\| \leq c h^{s-1} |u|_{s,\Omega} \]

Estimation \$L^2\$

Soient les hypothèse du théorème précédent, auxquelles nous ajoutons des hypothèses initial a des propriétés régularisantes Nous avons \$\lim_{h\rightarrow 0} \|u-u_h\|_{1,\Omega} = 0\$ et pour \$\frac{d}{2} < s < k+1, il existe une constante stem:[c\$ telle que \[ \forall h, \quad \|u-u_h\|{0,\Omega} \leq c h |u-u_h|{1,\Omega} \]

nous n’avons pas défini ce que sont ces propriétés régularisantes. Nous supposerons qu’elles sont vérifiées.
Exemple du Laplacien avec conditions de Dirichlet homogène en P1

\forall h, \quad \|u-u\|{0,\Omega} + h \|u-u_h\|{1,\Omega} \leq c h^2 \|f\|_{0,\Omega}

Convection dominated flows

Let \$\Omega \subset \mathbb{R}^d, d=1,2,3\$, consider this equation defined in \$\Omega\$, find \$u\$ a scalar field (e.g. temperature or concentration) such that

\$ \underbrace{\frac{\partial u}{\partial t}}_{time} \underbrace{- \nabla \cdot ( \epsilon \nabla u )}_{diffusion} + \underbrace{\mathbf{ \beta} \cdot \nabla u}_{advection/convection} + \underbrace{\mu u}_{reaction} = \underbrace{f}_{source},\quad \nabla \cdot \mathbf{ \beta} = 0\$

with the given data

  • \$\epsilon\$ is the diffusion coefficient (matrix \$d \times d\$), \$\epsilon(\mathbf{x}) > 0\$

  • \$\mathbf{\beta}\$ velocity and \$\nabla \cdot \mathbf{ \beta} \in L^2(\Omega)\$

  • \$\mu > 0\$ reaction (or absorption) coefficient.

From now on we consider only steady problems (\$\frac{\partial u}{\partial t}=0\$)

Applications

  • Transport of temperature

  • Transport of pollutants

  • Transport of chemical species (\$O_2, CO_2,...\$)

These equations are often coupled with fluid flow equations such as incompressible Navier-Stokes equations (air, blood, water…​) or darcy equations (porous media) to obtain the velocity field \$\mathbf{ \beta}\$. In this lesson we consider \$\mathbf{ \beta}\$ given.

Boundary conditions

Boundary conditions can be of 2 types on parts of the boundary or the whole boundary of the domain \$\Omega\$:

  • Dirichlet conditions : \$\label{eq:2} u|_{{\partial \Omega}_D} = g\$

  • Neumann conditions :

    • Outflow \$\label{eq:4} ( - \epsilon \nabla u + \mathbf{ \beta} u ) \cdot \mathbf{ n}|_{{\partial \Omega}_N} = (\mathbf{ \beta} u) \cdot \mathbf{ n}\$

    • (Heat) Flux \$\label{eq:5} ( - \epsilon \nabla u ) \cdot \mathbf{ n}|_{{\partial \Omega}_N} = Q \quad (e.g. \mathbf{\beta} = \mathbf{0})\$

Variational Formulation

The variational problem reads, Find \$u \in V\$ such that for all \$v\ in V\$

\$ a(u,v) \equiv \int_\Omega \epsilon \nabla u \cdot \nabla v + (\mathbf{ \beta} \cdot \nabla u ) v + \underbrace{\int_{\partial \Omega} (- \epsilon \nabla u ) \cdot \mathbf{ n}\ v}_{\mbox{apply boundary conditions}} = \int_\Omega f v \equiv \ell (v)\$

where

  • \$V\$ is an Hilbert space

  • we used the identity below for the integration by part

\$\nabla \cdot (\mathbf{ c} w ) = \mathbf{ c} \cdot \nabla w + \underbrace{w \nabla \cdot \mathbf{ c}}_{=0}\$
  • Suppose that \$-\frac{1}{2} \nabla \cdot ( \mathbf{ \beta} ) + \mu \geq 0\$ almost everywhere in \$\Omega\$ then we can show that \$a\$ is coercive provide \$\epsilon \geq \epsilon_0 > 0\$

\$a(v,v) \geq \alpha ||v||_{1,\Omega}, \quad \alpha = \frac{\epsilon_0}{1+C_\Omega^2}\$
  • \$a\$ is continuous, there exists \$M\$ a constant such that

\$|a(u,v)| \leq < M ||u||_{H^1{(\Omega)}} ||v||_{H^1{(\Omega)}}, \quad M = ||\mu||_{0,\Omega} + ||\epsilon||_{\infty,\Omega} + ||\mathbf{\beta}||_{\infty,\Omega}\$

Discrete formulation

Let \$V_h \subset V \equiv H^1(\Omega)\$ a discrete space, consider the standard Galerkin approximation on \$V_h\$ for . The problem reads

Find \$u_h \in V_h\$ such that for all \$v_h \in V_h\$ we have:

\[ \int_\Omega \epsilon \nabla u_h \cdot \nabla v_h + (\mathbf{ \beta} \cdot \nabla u_h ) v_h + \underbrace{\int_{\partial \Omega} (- \epsilon \nabla u_h ) \cdot \mathbf{ n}\ v_h}_{\mbox{apply boundary conditions}} = \int_\Omega f v_h \]

We can show that

\$||u_h||_{1,\Omega} \leq \frac{1}{\alpha} ||f||_{0,\Omega}, \quad ||\nabla u_h||_{1,\Omega} \leq \frac{C_\Omega}{\epsilon_0} ||f||_{0,\Omega},\$

which means that the Galerkin error inequality (best approximation) gives

\$||u-u_h||_V \leq \frac{M}{\alpha} \mathrm{inf}_{v_h \in V_h} ||u-v_h||_V\$

which means that given \$M\$ and \$\alpha\$, the estimate \$\epsilon_0 \rightarrow 0\$ In such case, the standard Galerkin method can yield to inaccurate solutions unless:

  • we use a very small \$h\$ (cost!!)

  • we use a stabilisation method (loss of \$\frac{1}{2}\$ in convergence rate)

Stabilisation methods for convection dominated flows

We now introduce

  • \$\mathrm{Pe}=\frac{|\mathbf{ \beta}|L}{2 \epsilon}\$ the global number. \$L\$ is the characteristic size of the domain.

  • the local Péclet number \$\mathrm{Pe}_h=\frac{|\mathbf{\beta}|h}{2 \epsilon}\$

The dominating convection and inacurrate behavior occurs when, on at least some cells, \$\mathrm{Pe} > 1\$ or \$\mathrm{Pe}_h > 1\$.

A remedy is to use a sufficiently small \$h\$ same to ensure \$\mathrm{Pe}_h < 1\$. For example if \$|\mathbf{beta}| = 1\$ and \$\epsilon=5e-4\$, we should take \$h=1e-4\$.
Another remedy is to use a different approximation space for the unknown \$u_h\$ and the test functions \$v_h\$. We talk about Petrov-Galerkin approximation.

Find \$u_h \in V_h\$ such that

\[ a_h(u_h,v_h) = \ell_h(v_h),\quad \forall v_h \in V_h \] with

\[ a_h(u_h,v_h) = a(u_h,v_h) + b_h(u_h,v_h),\quad \ell_h(v_h) = \ell(v_h) + g_h(v_h),\quad \].

The purpose of \$b_h\$ and \$g_h\$ is to eliminate(or reduce) the numerical oscillations produced by the standard Galerkin method: they are called stabilisation methods and depend on \$h\$. The term stabilisation method is not exact. Indeed the Galerkin method is already stable (i.e. continuity). Here it is to be understood as the aim of reducing (or elimination) numerical oscillations when \$\mathrm{Pe} > 1\$.

Without doing anything wiggles occur.

There are remedies so called stabilisation techniques, here some some examples:

  • Artificial diffusion (streamline diffusion) (SDFEM)

  • Galerkin Least Squares method (GaLS)

  • Streamline Upwind Petrov Galerkin (SUPG)

  • Continuous Interior Penalty methods (CIP)

Artificial diffusion (or streamline diffusion) (SDFEM)

Method The method consists in adding an

\$\epsilon_h =\epsilon(1+\phi(\mathrm{Pe}))\$

with \$\phi(\mathrm{Pe}) \rightarrow 0\$ as \$h \rightarrow 0\$, e.g. \$\phi(\mathrm{Pe}) = \mathrm{Pe}-1+B(2*\mathrm{Pe})\$ where \$B\$ is the so-called Bernoulli function \$B(t) = \frac{t}{e^t-1}\$ if \$t > 0\$ and \$B(0) = 1\$ (also exponential fitting scheme)

\$ b_h(u_h,v_h) = \int_\Omega \epsilon \Phi(\mathrm{Pe}) \nabla u_h \cdot \nabla v_h, \quad g_h(v_h) = 0\$

for a given \$\epsilon\$ and for \$h\$ tending to \$0\$, we have for \$u \in H^{r+1}(\Omega)\$

\[ ||u-u_h||{1,\Omega} \leq C_1 \Big[ h^r||u||{r+1,\Omega} + \phi(\mathrm{Pe})||u||_{1,\Omega}\Big] \] and for a given \$h\$ and \$\epsilon\$ tending to 0,

\[ ||u-u_h||{1,\Omega} \leq C_1 \Big[ h^{r-1}||u||{r+1,\Omega} + ||u||_{1,\Omega}\Big] \] If \$\phi(\mathrm{Pe})=\frac{|\mathbf{ \beta}|h}{2 \epsilon}\$, the convergence is linear, with the exponential fitting scheme it is quadratic if \$r \geq 2\$.

GaLS and SUPG

First we decompose our operators into a symmetric (\$<Lu,v> = <u,Lv>\$ and skew symmetric (\$<L u, v> = -<u,L v>\$) contributions, we start with

\$ L u = -\epsilon \Delta u + \nabla \cdot (\mathbf{ \beta} u ) + \mu u\$
\$L u = \underbrace{-\epsilon \Delta u + \Big[ \mu + \frac{1}{2} \nabla \cdot \mathbf{ \beta} \Big] u}_{L_S u} + \underbrace{\frac{1}{2}\Big[ \nabla \cdot ( \mathbf{ \beta} u) + \mathbf{ \beta} \cdot \nabla u \Big]}_{L_{SS} u}\$
Consistent schemes

We say that a method is consistent when adding a term to a problem such as:

Find \$u_h \in V_h\$ such that

\[ a(u_h,v_h) + \mathcal{L}_h(u_h,f;v_h) = (f,v_h), \quad \forall v_h \in V_h\] the term added statisfies

\[ \mathcal{L}_h(u,f;v_h) = 0, \forall v_h \in V_h \]

Choice for consistent methods

A possible choice for \$\mathcal{L}_h\$ is the following

\$ \mathcal{L}_h(u_h,f;v_h) = \mathcal{L}^{(\rho)}_h(u_h,f;v_h) = \sum_{K \in \mathcal{T}_h} \delta (L u_h - f, \mathcal{S}^{(\rho)}_K(v_h))_{0,\Omega}\$

where

  • \$(\cdot,\cdot)_{0,\Omega}\$ is the \$L^2\$ scalar product

  • \$\rho\$ and \$\delta\$ are parameters

and we have set

\$\mathcal{S}^{(\rho)}_K(v_h) = \frac{h_K}{|\mathbf{\beta}|}\Big[ L_{SS} v_h + \rho L_S v_h\Big]\$
Galerkin Least-Square

if \$\rho = 1\$ we have the Galerkin Least Square method (GaLS)

\$\mathcal{S}^{(\rho)}_K(v_h) = \frac{h_K}{|\mathbf{ \beta}|}\Big[ L v_h\Big]\$
Streamline Upwind Petrov-Galerkin

if \$\rho = 0\$ we have the Streamline Upwind Petrov-Galerkin (SUPG)

\$\mathcal{S}^{(0)}_K(v_h) = \frac{h_K}{|\mathbf{ \beta}|}\Big[ L_{SS} v_h\Big]\$
Douglas and Wang

if \$\rho = -1\$ we have the Douglas and Wang (DW)

\$\mathcal{S}^{(-1)}_K(v_h) = \frac{h_K}{|\mathbf{ \beta}|}\Big[ (L_{SS} -L_S )v_h\Big]\$

We define the \$\rho\$ Norm

\$||v||_{(\rho)} = \Big\{\epsilon ||\nabla u||^2_{0,\Omega} + ||\sqrt{\gamma} v||^2_{0,\Omega} + \sum_{K \in \mathcal{T_h}} \delta \Big( (L_{SS}+\rho L_S )v, \mathcal{S}^{(\rho)}_K(v) \Big)_{0,\Omega} \Big\}^{1/2}\$

where \$\gamma\$ is a positive constant such that \$-\frac{1}{2} \nabla \cdot \mathbf{\beta} + \mu \geq \gamma > 0\$

We have the following result

if \$u \in H^{k+1}(\Omega)\$, then the following error estimates hold:

\[ {\|u-u_h\|{(\rho)}} \leq C {h^{k+1/2}} |u|{k+1,\Omega} \]

GaLS

In practice for GaLS (\$\rho = 1\$) we take \$\delta\$ such that

\$\delta(h_K,\epsilon) \frac{h_K}{|\mathbf{ \beta}|} = \Big( \frac{1}{h_K} + \frac{\epsilon}{h^2_K} \Big)^{-1}\$

and we can prove the following estimates if \$u\in H^{k+1}(\Omega)\$,

\$\forall \epsilon \quad {\|u-u_h\|_{0,\Omega}} \leq c {h^{k+1/2}} \|u\|_{k+1,\Omega}\$
\$\forall \epsilon \geq c h \quad {\|u-u_h\|_{1,\Omega}} \leq c {h^{k}} \|u\|_{k+1,\Omega}\$

and finally if the family \$\{\mathcal{T}_h\}_{h > 0}\$ is quasi-uniform and \$\epsilon \leq c h \$, then

\$\| \beta \cdot \nabla (u -u_h) \|_{0,\Omega} \leq c h^k \|u \|_{k+1,\Omega}\$

Continuous Interior Penalty

In the continuous interior penalty we add the following term

\$\sum_{F \in \Gamma_\mathrm{int} } \int_{F} \eta\ h_F^2\ |\mathbf{ \beta} \cdot \mathbf{n}|\ \jump{\nabla u} \jump{\nabla v}\$

where

  • \$\Gamma_\mathrm{int}\$ is the set of internal faces

  • the \$\mathrm{Pe}>>1\$ (typically it is applied to all internal faces)

  • we have

\jump{\nabla u} = \nabla u \cdot \mathbf{n}|_1 + \nabla u \cdot \mathbf{n}|_2

is the so called jump of \$\nabla u\$(scalar valued) across the face.

In the case of scalar valued functions

\$ \jump{u} = u \mathbf{n}|_1 + u \mathbf{n}|_2\$

Choice for \$\eta\$ \$\eta\$ can be taken in the range \$[1e-2;1e-1\$]. A typical value is \$\eta=2.5e-2\$. A similar error estimate \$O(h^{r+1/2})\$ holds for CIP.

Example CIP

// define the stabilisation coefficient expression
auto stab_coeff = (eta*abs(idv(beta))*abs(trans(N())*idv(beta)))*vf::pow(hFace(),2.0));

// assemble the stabilisation operator
form2( Xh, Xh, M ) +=
 integrate( internalfaces(Xh->mesh()), // faces of the mesh
            stab_coeff*(trans(jumpt(gradt(u)))*jump(grad(v))));

Élasticité Linéaire

On s’intéresse dans cette section à l’approximation par éléments finis de problèmes de mécanique des milieux continus en 3D.

Soit \$\mathbf{f}: \Omega \rightarrow \mathbb{R}^3\$ la charge extérieur s’appliquant au domaine \$\Omega\$.

On note \$\disp: \Omega \rightarrow \mathbb{R}^3\$ le déplacement de la structure induit par cette charge \$\mathbf{f}\$.

En supposant que les déformations soient petites pour être modélisées dans le cadre de l’elasticité linéaire, on a la relation suivante à l’équilibre:

\[ \nabla \cdot \stresst + \mathbf{f} = \mathbf{0} \mbox{ dans } \Omega \] où \$\stresst: \Omega \rightarrow \RR^{d\times d}\$ est le tenseur des contraintes défini par la relation

\[ \stresst = \lambda \tr(\deformt) \Id + 2 \mu \deformt \] et \$\deformt : \Omega \rightarrow \RR^{d\times d}\$ est le tenseur des déformations défini par

\[ \deformt = \frac{1}{2} \left( \nabla \disp + \nabla \disp^T \right), \] \$\lambda\$ et \$\mu\$ les coefficients de Lamé et \$\Id\$ la matrice identité de \$\RR^{d\times d}\$.

On a alors

\[ \stresst = \lambda( \nabla \cdot \disp ) \Id + \mu( \nabla \disp + \nabla \disp^T) \]

Coefficients de Lamé

Ils sont des coefficients phénoménologiques contraints par les relations suivants:

  • \$\mu >0\$

  • \$\lambda + \frac{2}{3} \mu \ge 0\$

Dans ce qui suit, on supposera que \$\lambda \ge 0\$ et ces coefficients constants.

D’un point de vue pratique, ces coefficients sont obtenus par les module d’Young \$E\$ et coefficient de Poisson \$\nu\$ tels que

\[ \lambda = \frac{E \nu}{( 1+\nu )*( 1-2 \nu )} , \quad \mu =\frac{E}{2 ( 1+\nu )} \]