Equation de diffusion

Christophe Prud'homme <christophe.prudhomme@cemosis.fr>, Laurent Navoret

Problème de diffusion non linéaire

Equation de Poisson non-linéaire

\[-\nabla \cdot(q(u) \nabla u)=f\]

Le coefficient \(q(u)\) rend l’équation non linéaire (sauf si \(q(u)\) est constant en \(u\) ).

Verification

Pour pouvoir vérifier facilement notre implémentation, nous choisissons le domaine, \(q(u), f,\) et les conditions aux limites de telle sorte que nous ayons une solution simple et exacte \(u\).

Soit \(\Omega\) l’hypercube unitaire \([0,1]^{d}\) en \(d\) dimensions, \(q(u)=(1+u)^{m}, f=0, u=0\) pour \(x_{0}=0, u=1\) pour \(x_{0}=1,\) et \(\partial u / \partial n=0\) à toutes les autres frontières \(x_{i}=0\) et \(x_{i}=1, i=1, \ldots, d-1\).

Les coordonnées sont maintenant représentées par les symboles \(x_{0}, \ldots, x_{d-1} .\)

La solution exacte est alors

\[u\left(x_{0}, \ldots, x_{d}\right)=\left(\left(2^{m+1}-1\right) x_{0}+1\right)^{1 /(m+1)}-1\]

Formulation variationelle

La formulation variationnelle de notre problème modèle est la suivante : trouver \(u \in V\) tel que

\[F(u ; v)=0 \quad \forall v \in \hat{V}\]

où

\[F(u ; v)=\int_{\Omega} q(u) \nabla u \cdot \nabla v \mathrm{~d} x\]

\[\begin{array}{l} \hat{V}=\left\{v \in H^{1}(\Omega) : v=0 \text{ sur } x_{0}=0 \text{ et } x_{0}=1\right\} \\ V=\left\{v \in H^{1}(\Omega) : v=0 \text{ sur } x_{0}=0 \text{ et } v=1 \text{ sur } x_{0}=1\right\} \end{array}\]

Discrete Problem

Le problème discret se pose comme d’habitude en restreignant \(V\) et \(\hat{V}\) à une paire d’espaces discrets. Comme souvent, nous omettons tout indice sur les espaces discrets et disons simplement que \(V\) et \(\hat{V}\) sont choisis en dimension finie selon un certain type de maille et d’élément. Le problème non linéaire est alors le suivant : trouver \(u \in V\) tel que

\[F(u ; v)=0 \quad \forall v \in \hat{V}\]

avec \(u=\sum_{j=1}^{N} U_{j} \phi_{j}\). Puisque \(F\) est une fonction non linéaire de \(u\) donne lieu à un système d’équations algébriques non linéaires.

Stratégies de résolution

3 Stratégies

Picard
Newton au niveau algébrique
Newton au niveau EDP

Picard

L'itération de Picard est un moyen facile de traiter les EDP non linéaires : nous utilisons simplement une solution antérieure connue dans les termes non linéaires de sorte que ces termes deviennent linéaires dans l’inconnue \(u\).

Cette stratégie est également connue comme la méthode des substitutions successives. Pour notre problème particulier, nous utilisons une solution antérieure connue dans le coefficient \(q(u)\). Plus précisément, étant donné une solution \(u^{k}\) de l’itération \(k\), nous cherchons une nouvelle solution (que nous espérons améliorée) \(u^{k+1}\) dans l’itération \(k+1\) telle que \(u^{k+1}\) résout le problème linéaire

\[\nabla \cdot\left(q\left(u^{k}\right) \nabla u^{k+1}\right)=0, \quad k=0,1, \ldots\]

Les itérations nécessitent une estimation initiale \(u^{0}.\) On espère que \(u^{k} \rightarrow u\) comme \(k \rightarrow \infty,\) et que \(u^{k+1}\) soit suffisamment proche de la solution exacte \(u\) du problème discret après seulement quelques itérations.

Picard (2)

Nous pouvons facilement formuler un problème variationnel pour \(u^{k+1}\). De manière équivalente, nous pouvons approximer \(q(u)\) par \(q\left(u^{k}\right)\) pour obtenir le même problème variationnel linéaire. Dans les deux cas, le problème consiste à chercher \(u^{k+1} \in V\) tel que

\[\tilde{F}\left(u^{k+1} ; v\right)=0 \quad \forall v \in \hat{V}, \quad k=0,1, \ldots\]

\[\mbox{avec } \tilde{F}\left(u^{k+1} ; v\right)=\int_{\Omega} q\left(u^{k}\right) \nabla u^{k+1} \cdot \nabla v \mathrm{~d} x\]

Puisqu’il s’agit d’un problème linéaire dans l’inconnue \(u^{k+1}\), ou de manière équivalente

\[a\left(u^{k+1}, v\right)=L(v)\\ \mbox{avec } a(u, v) =\int_{\Omega} q\left(u^{k}\right) \nabla u \cdot \nabla v \mathrm{~d} x \quad L(v) =0\]

Picard (Iterations)

Les itérations peuvent être arrêtées lorsque

\(\epsilon \equiv \| u^{k+1}-u^{k}||<\texttt{tol}\), où \(\texttt{tol}\) est petit, disons \(10^{-5}\),
le nombre d’itérations dépasse une certaine limite critique.

Dans ce dernier cas, on constatera une divergence de la méthode ou une convergence lente inacceptable.

Newton algébrique 1

Après avoir discrétisé notre problème d’EDP non linéaire, nous pouvons utiliser la méthode de Newton pour résoudre le système d’équations algébriques non linéaires.

A partir du problème variationnel continu, la version discrète résulte en un système d’équations pour les paramètres inconnus \(U_{1}, \ldots, U_{N}\) (en insérant \(u=\sum_{j=1}^{N} U_{j} \phi_{j}\) et \(v=\hat{\phi}_{i}\) dans [nl-discrete]):

\[F_{i}\left(U_{1}, \ldots, U_{N}\right) \equiv \sum_{j=1}^{N} \int_{\Omega}\left(q\left(\sum_{\ell=1}^{N} U_{\ell} \phi_{\ell}\right) \nabla \phi_{j} U_{j}\right) \cdot \nabla \hat{\phi}_{i} \mathrm{~d} x=0, \quad i=1, \ldots, N\]

Newton algébrique 2

La méthode de Newton pour le système \(F_{i}\left(U_{1}, \ldots, U_{j}\right)=0, i=1, \ldots, N\) peut être formulée comme suit

\[\begin{aligned} \sum_{j=1}^{N} \frac{\partial}{\partial U_{j}} F_{i}\left(U_{1}^{k}, \ldots, U_{N}^{k}\right) \delta U_{j} &=-F_{i}\left(U_{1}^{k}, \ldots, U_{N}^{k}\right), \quad i=1, \ldots, N \\ U_{j}^{k+1} &=U_{j}^{k}+\omega \delta U_{j}, \quad j=1, \ldots, N \end{aligned}\]

où \(\omega \in[0,1\)] est un paramètre de relaxation, et \(k\) est un indice d’itération.

Remarques

Une estimation initiale \(u^{0}\) doit être fournie pour démarrer l’algorithme.
La méthode originale de Newton a \(\omega=1\), mais dans les problèmes où il est difficile d’obtenir la convergence, ce qu’on appelle la sous-relaxation avec \(\omega<1\) peut aider.

Newton algébrique 3

Nous avons besoin, dans un programme, de calculer la matrice jacobienne \(\partial F_{i} / \partial U_{j}\) et le vecteur du côté droit \(-F_{i}\). La dérivée \(\partial F_{i} / \partial U_{j}\) s’écrit

\[\int_{\Omega}\left[q^{\prime}\left(\sum_{\ell=1}^{N} U_{\ell}^{k} \phi_{\ell}\right) \phi_{j} \nabla\left(\sum_{j=1}^{N} U_{j}^{k} \phi_{j}\right) \cdot \nabla \hat{\phi}_{i}+q\left(\sum_{\ell=1}^{N} U_{\ell}^{k} \phi_{\ell}\right) \nabla \phi_{j} \cdot \nabla \hat{\phi}_{i}\right] \mathrm{d} x\]

Les résultats suivants ont été utilisés pour obtenir l’équation ci-dessus:

\[\frac{\partial u}{\partial U_{j}}=\frac{\partial}{\partial U_{j}} \sum_{j=1}^{N} U_{j} \phi_{j}=\phi_{j}, \quad \frac{\partial}{\partial U_{j}} \nabla u=\nabla \phi_{j}, \quad \frac{\partial}{\partial U_{j}} q(u)=q^{\prime}(u) \phi_{j}\]

Newton algébrique 4

Nous pouvons reformuler la matrice jacobienne en introduisant la notation courte \(u^{k}=\sum_{j=1}^{N} U_{j}^{k} \phi_{j}\) ;

\[\frac{\partial F_{i}}{\partial U_{j}}=\int_{\Omega}\left[q^{\prime}\left(u^{k}\right) \phi_{j} \nabla u^{k} \cdot \nabla \hat{\phi}_{i}+q\left(u^{k}\right) \nabla \phi_{j} \cdot \nabla \hat{\phi}_{i}\right] \mathrm{d} x\]

Afin de faire calculer cette matrice, nous devons formuler un problème variationnel correspondant. En regardant le système d’équations linéaires de la méthode de Newton,

\[\sum_{j=1}^{N} \frac{\partial F_{i}}{\partial U_{j}} \delta U_{j}=-F_{i}, \quad i=1, \ldots, N\]

nous pouvons introduire \(v\) comme fonction de test générale remplaçant \(\hat{\phi}_{i}\), et nous pouvons identifier l’inconnue \(\delta u=\sum_{j=1}^{N} \delta U_{j} \phi_{j}\).

Newton algébrique 5

À partir du système linéaire, nous pouvons maintenant "revenir en arrière" pour construire la forme faible discrète correspondante.

\[\int_{\Omega}\left[q^{\prime}\left(u^{k}\right) \delta u \nabla u^{k} \cdot \nabla v+q\left(u^{k}\right) \nabla \delta u \cdot \nabla v\right] \mathrm{d} x=-\int_{\Omega} q\left(u^{k}\right) \nabla u^{k} \cdot \nabla v \mathrm{~d} x\]

L’équation s’adapte à la forme standard \(a(\delta u, v)=L(v)\) avec

\[\begin{aligned} a(\delta u, v) =\int_{\Omega}\left[q^{\prime}\left(u^{k}\right) \delta u \nabla u^{k} \cdot \nabla v+q\left(u^{k}\right) \nabla \delta u \cdot \nabla v\right] \mathrm{d} x, \quad L(v) =-\int_{\Omega} q\left(u^{k}\right) \nabla u^{k} \cdot \nabla v \mathrm{~d} x \end{aligned}\]

Remarques

La solution précédente \(u^{k}\) remplace \(u\) dans les formules lors du calcul de la matrice \(\partial F_{i} / \partial U_{j}\) et le vecteur \(F_{i}\) pour le système linéaire à chaque itération de Newton.
Afin d’obtenir une bonne estimation initiale \(u^{0}\), nous pouvons résoudre un problème linéaire simplifié, typiquement avec \(q(u)=1,\) qui donne l’équation de Laplace standard \(\Delta u^{0}=0.\)

Newton EDP (1)

Bien que la méthode de Newton dans les problèmes d’EDP soit normalement formulée au niveau de l’algèbre linéaire, c’est-à-dire comme une méthode de résolution de systèmes d’équations algébriques non linéaires, nous pouvons également formuler la méthode au niveau des EDP.

Cette approche permet de linéariser les EDP avant de les discrétiser.

Étant donné une approximation du champ de solution, \(u^{k}\), nous recherchons une perturbation \(\delta u\) telle que

\[u^{k+1}=u^{k}+\delta u\]

satisfasse l’EDP non linéaire. Cependant, le problème pour \(\delta u\) est toujours non linéaire et on ne gagne rien.

Newton EDP (2)

L’idée est donc de supposer que \(\delta u\) est suffisamment petit pour que l’on puisse linéariser le problème par rapport à \(\delta u\).

En insérant \(u^{k+1}\) dans l’EDP, on linéarise le terme \(q\) comme suit

\[q\left(u^{k+1}\right)=q\left(u^{k}\right)+q^{\prime}\left(u^{k}\right) \delta u+\mathcal{O}\left((\delta u)^{2}\right) \approx q\left(u^{k}\right)+q^{\prime}\left(u^{k}\right) \delta u\]

et en abandonnant les autres termes non linéaires dans \(\delta u,\) on obtient

\[\nabla \cdot\left(q\left(u^{k}\right) \nabla u^{k}\right)+\nabla \cdot\left(q\left(u^{k}\right) \nabla \delta u\right)+\nabla \cdot\left(q^{\prime}\left(u^{k}\right) \delta u \nabla u^{k}\right)=0\]

Nous pouvons rassembler les termes avec l’inconnue \(\delta u\) sur le côté gauche,

\[\nabla \cdot\left(q\left(u^{k}\right) \nabla \delta u\right)+\nabla \cdot\left(q^{\prime}\left(u^{k}\right) \delta u \nabla u^{k}\right)=-\nabla \cdot\left(q\left(u^{k}\right) \nabla u^{k}\right)\]

Newton EDP (3)

La forme faible de cette EDP est dérivée en la multipliant par une fonction test \(v\) et en l’intégrant sur \(\Omega\), en intégrant les dérivées du second ordre par parties :

\[\int_{\Omega}\left(q\left(u^{k}\right) \nabla \delta u \cdot \nabla v+q^{\prime}\left(u^{k}\right) \delta u \nabla u^{k} \cdot \nabla v\right) \mathrm{d} x=-\int_{\Omega} q\left(u^{k}\right) \nabla u^{k} \cdot \nabla v \mathrm{~d} x\]

Le problème variationnel est le suivant : trouver \(\delta u \in V\) tel que \(a(\delta u, v)=L(v)\) pour tout \(v \in \hat{V},\) où

\[\begin{aligned} a(\delta u, v) &=\int_{\Omega}\left(q\left(u^{k}\right) \nabla \delta u \cdot \nabla v+q^{\prime}\left(u^{k}\right) \delta u \nabla u^{k} \cdot \nabla v\right) \mathrm{d} x \\ L(v) &=-\int_{\Omega} q\left(u^{k}\right) \nabla u^{k} \cdot \nabla v \mathrm{~d} x \end{aligned}\]

Les espaces de fonctions \(V\) et \(\hat{V}\), étant continus ou discrets, sont comme dans le problème de Poisson linéaire.

Newton EDP (4)

Nous devons fournir une estimation initiale, par exemple, la solution de l’EDP avec \(q(u)=1\).

La forme faible correspondante \(a_{0}\left(u^{0}, v\right)=L_{0}(v)\) possède

\[a_{0}(u, v)=\int_{\Omega} \nabla u \cdot \nabla v \mathrm{~d} x, \quad L(v)=0\]

Ensuite, nous entrons dans une boucle et résolvons \(a(\delta u, v)=L(v)\) pour \(\delta u\) et calculons une nouvelle approximation \(u^{k+1}=u^{k}+\delta u\).

Remarques

\(\delta u\) est une correction, donc si \(u^{0}\) satisfait les conditions de Dirichlet prescrites sur une partie \(\Gamma_{\mathrm{D}}\) de la frontière, nous devons exiger \(\delta u=0\) sur \(\Gamma_{\mathrm{D}}\).
Le critère d’arrêt de la boucle se fait sur l’incrément \(\|\delta u\|\)
L’algorithme (et le problème linéaire à résoudre) est le même que dans le cas algébrique

Equation de diffusion dépendant du temps

\[\begin{split} \partial_{t} u-\Delta u&=f, \quad \mbox{ sur } \Omega \times \mathbb{R}^{+}\\ u(x,t=0)&=u_{0}\\ u(x,t)&=g(x,t), \quad \mbox{ sur } \partial \Omega_{D} \times \mathbb{R}^{+} \end{split}\]

\(u: \Omega \times \mathbb{R}^{+} \rightarrow \mathbb{R}:\) densité
\(u_{°}: \Omega \rightarrow \mathbb{R}:\) donnée initiale
\(f: \Omega \times \mathbb{R}^{+} \rightarrow \mathbb{R}:\) terme source
\(g: \partial \Omega \rightarrow \mathbb{R}:\) valeur au bord

Formulation

\(\rightarrow\) problème parabolique \(\rightarrow\) Bilan de conservation

\[\frac{d}{d t} \int_{\Omega} u=\int_{\partial \Omega}(n \cdot \nabla u)+\int_{\Omega} f\]

\(\rightarrow\) Dissipation de la norme \(L^{2}\)

\[\frac{d}{d t} \int_{\Omega} u^{2}=-\int_{\Omega}|\nabla u|^{2}+\int_{\partial_{\Omega}} u \partial_{n} u+\int_{\Omega} f u\]

Approche

Une approche simple pour résoudre les EDP dépendant du temps par La méthode des éléments finis consiste à discrétiser d’abord la dérivée en temps par les différences finies. Cela donne un ensemble récursif de problèmes stationnaires, puis on transforme chaque problème stationnaire en un problème variationnel.

Soit l’exposant \(k\) qui désigne une quantité au temps \(t^k\) où \(k\) est un nombre entier comptant les niveaux de temps. Par exemple, \(u^k\) signifie \(u\) au niveau du temps \(k\).

Semi-discrétisation en temps

Discrétisation en temps: \(\Delta t>0, t_{n}=n \Delta t .\) On note \(u^{n}(x):=u\left(t_{n}, x\right)\) Schéma d’Euler (implicite):

\[\frac{u^{n}-u^{n-1}}{\Delta t}=\Delta u^{n}+f^{n}\]

\(\rightarrow\) "method of lines" \(\rightarrow\) pour trouver \(u^{n},\) on résout un problème elliptique (stationnaire)

\[u^{n}-\Delta t \Delta u^{n}=u^{n-1}+\Delta t f^{n}\]

Résolution en espace

Problème au temps \(t^{n}\)

\[u^{n}-\Delta t \Delta u^{n}=u^{n-1}+\Delta t f^{n}\]

\(V=\left\{v \in H^{1}(\Omega), v=g^{n} \text { sur } \partial \Omega\right\}\) \(V_{0}=\left\{v \in H^{1}(\Omega), v=0 \text { sur } \partial \Omega\right\}\) Problème variationnel

\[\text { Trouver } u^{n} \in V, \quad a\left(u^{n}, v\right)=\ell(v), \quad \forall v \in V_{0}, \mbox{ avec }\]

\[a(u, v)=\int_{\Omega} u v+\Delta t \nabla u \cdot \nabla v,\quad \ell(v)=\int_{\Omega}\left(u^{n-1}+\Delta t f^{n}\right) v\]

\(\rightarrow\) Méthode d’éléments finis

Remarques

Pour le calcul de \(u^{0}\) :

interpolation: \(u^{0}=\sum_{j=1}^{N} u_{0}\left(x_{i}\right) \varphi_{i}\)
projection \(L^{2}: u^{0}\) est solution du problème variationnel

\[\int_{\Omega} u v=\int_{\Omega} u_{0} v, \quad \forall v \in V_{h}\]

Autre schéma numérique pour la discrêtisation de la dérivée en temps : Crank-Nicolson (ordre 2 en temps)

Advection Diffusion Reaction

Problèmes type

Problème de convection-diffusion

\[\begin{array}{rlr} \partial_{t} u+\nabla \cdot(u \beta)-\varepsilon \Delta u+\mu u & =f & \text { sur } \mathbb{R}^{+} \times \Omega\\ u & =0 & \text { sur } \mathbb{R}^{+} \times \partial \Omega \end{array}\]

\(u: \mathbb{R}^{+} \times \Omega \rightarrow \mathbb{R}^{d}:\) température ou concentration
\(\beta: \mathbb{R}^{+} \times \Omega \rightarrow \mathbb{R}^{d}:\) vitesse d’advection

Problèmes type

Problème de Navier-Stokes

\[\begin{array}{rlr} \partial_{t} u+u \cdot \nabla u-\nu \Delta u+\nabla p & =f & & \text { sur } \mathbb{R}^{+} \times \Omega \\ \nabla \cdot u & =0 & & \text { sur } \mathbb{R}^{+} \times \Omega \\ u & =0 & & \text { sur } \mathbb{R}^{+} \times \partial \Omega \end{array}\]

\(u: \mathbb{R}^{+} \times \Omega \rightarrow \mathbb{R}^{d}:\) champ de vitesse
\(p: \mathbb{R}^{+} \times \Omega \rightarrow \mathbb{R}:\) pression

advection non-linéaire

Advection diffusion stationnaire

\[\begin{aligned} \nabla \cdot(u \beta)-\varepsilon \Delta u+\mu u &=f & & \text { sur } \mathbb{R}^{+} \times \Omega \\ u &=0 & & \text { sur } \mathbb{R}^{+} \times \partial \Omega \end{aligned}\]

Proposition

\[u \in H^{1}(\Omega) \text { et } \beta \in H^{1}(\Omega)^{d}\quad \nabla \cdot(u \beta)=u \nabla \cdot \beta+\beta \cdot \nabla u\]

1) Formulation variationnelle

\[\text { Trouver } u \in V, \quad a(u, v)=\ell(v), \quad \forall v \in V\equiv H_{0}^{1}(\Omega)\\ \begin{array}{l} a(u, v)=\int_{\Omega} \varepsilon \nabla u \cdot \nabla v-\int_{\Omega} u(\beta \cdot \nabla v)+\int_{\Omega} \mu u v \\ \ell(v)=\int_{\Omega} f v \end{array}\]

Résolution

Proposition: \(f \in L^{2}, \nabla \cdot \beta \in L^{2}\) et \(\beta \in L^{\infty}, \varepsilon, \mu>0 .\) On suppose de plus que \(\mu+\frac{1}{2} \nabla \cdot \beta \geqslant 0\) alors \(a\) est coercif : \(a(v, v) \geqslant \kappa\|v\|_{H^{1}}^{2},\) avec \(\kappa=\varepsilon\) \(a\) est continu \(: a(u, v) \leqslant M\|u\|_{H^{1}}\|v\|_{H^{1}},\) avec \(M=\left(\mu+\varepsilon+\|\beta\|_{\infty}\right)\) On en déduit qu’il existe une unique solution u au problème variationnel.
Definition ( Problème variationnel discret)

\[\begin{array}{c} \text { Soit } V_{h} \subset V \text{, Trouver } u_{h} \in V_{h}, \quad a\left(u_{h}, v_{h}\right)=\ell\left(v_{h}\right), \quad \forall v_{h} \in V_{h} \\ a\left(u_{h}, v_{h}\right)=\int_{\Omega} \varepsilon \nabla u_{h} \cdot \nabla v_{h}-\int_{\Omega} u_{h}\left(\beta \cdot \nabla v_{h}\right)+\int_{\Omega} \mu u_{h} v_{h} \\ \ell\left(v_{h}\right)=\int_{\Omega} f v_{h} \end{array}\]

2) Choix du maillage
3) Choix des éléments fins \(\left(P_{k}\right)\)

Problème

Convergence: Preuve de la convergence :

\[\left\|u-u_{h}\right\| \leqslant \frac{M}{\kappa} \operatorname{dist}\left(u, V_{h}\right)=\frac{\left(\mu+\varepsilon+\|\beta\|_{\infty}\right)}{C \varepsilon} \operatorname{dist}\left(u, V_{h}\right)\]

On a ensuite :

\[\left\|u_{h}\right\| \leqslant \frac{1}{\varepsilon}\|f\|_{L^{2}}, \quad\left\|\nabla u_{h}\right\| \leqslant \frac{1}{\varepsilon}\|f\|_{L^{2}}\]

Limite faible diffusion (convection dominante): \(\varepsilon \rightarrow 0\) \(\rightarrow\) solution et gradient très grands

Nombre de Péclet

Definition

Nombre de Péclet: \(\mathrm{Pe}=\frac{\|\beta\|_{\infty} L}{2 \varepsilon}(L=\) longueur caractéristique du domaine)
Nombre de Péclet local à une maille: \(\mathrm{Pe}_{K}=\frac{\|\beta\|_{\infty} h_{K}}{2 \varepsilon}\)

Pour limiter les oscillations

réduire l’erreur en prenant un très petit paramètre de discrétisation \(h\)
méthodes de stabilisation :
- Diffusion artificielle
- Galerkin Least-Square (GLS)
- Streamline Upwind Petrov Galerkin (SUPG)
- Continuous Interior Penalty methods (CIP)
- ..

Méthode 1: Galerkin Least Square (GaLS)

Problème : \(L u=f\) avec

\[L u=\nabla \cdot(u \beta)-\varepsilon \Delta u+\mu u\]

Formulation moindre carré :

\[\text { Trouver } u \in V, \quad(L u, L v)=(f, L v), \quad \forall v \in V\]

Cela revient à résoudre \(L^{T} L u=L^{T} f .\) L’unique solution de ce problème est le minimum de la fonction

\[J(u)=\frac{1}{2}(L u, L u)-(f, L u)=\frac{1}{2}\|L u-f\|^{2}-\frac{1}{2}\|f\|^{2}\]

Problème mal posé : pour \(v \in H^{1}\), \(L v \notin L^{2}(\Omega)\)

Galerkin Least Square (GaLS)

mais pour tout \(v_{h} \in V_{h}\) et tout \(K \in \mathcal{T}_{h}, L v_{h \mid K} \in L^{2}\).

Formulation variationnelle GLS: Trouver \(u \in V_{h},\)

\[a(u, v)+\underbrace{\sum_{K \in \mathcal{T}_{h}}\left(L u_{\mid K}, \delta_{K} L v_{\mid K}\right)}_{a_{h}(u, v)}=\ell(v)+\underbrace{\sum_{K \in \mathcal{T}_{h}}\left(f, \delta_{K} L v_{\mid K}\right)}_{\ell_{h}(v)}, \quad \forall v \in V_{h}\]

Remarques

\(\rightarrow\) combinaison entre méthode de Galerkin standard et méthode de moindre carré
\(\rightarrow\) méthode fortement consistante : la solution \(u\) vérifie la formulation variationnelle
\(\rightarrow \delta_{K}:\) paramètre tel que \(\delta_{K}=O\left(h_{K} /\|\beta\|\right)\) pour \(\mathrm{Pe}_{K} \geqslant 1\)
\(\delta_{K}=O\left(h_{K}^{2} / \varepsilon\right)\) pour \(\mathrm{Pe}_{K} \leqslant 1\)
\(\rightarrow a_{h}\) et \(\ell_{h}\) dépendent de \(\left(h_{K}\right)\)

Convergence

Norme

\[\|u\|_{\mathrm{GLS}}^{2}=\varepsilon\|\nabla u\|_{L^{2}}^{2}+\|\sqrt{\tilde{\mu}} u\|_{L^{2}}^{2}+\sum_{K \in \mathcal{T}_{h}}\left(L u_{\mid K}, \delta_{K} L u_{\mid K}\right), \quad \tilde{\mu}=\mu+\frac{1}{2} \nabla \cdot \beta \geqslant \mu_{0}>0\]

Proposition: Soit \(f \in L^{2}(\Omega)\) et \(u \in V\) la solution du problème. Soit \(u_{h} \in V_{h}\) la solution approchée du problème d’advection-diffusion. Si \(u \in H^{k+1}(\Omega),\) alors avec \(\delta_{K}=\delta h_{K} /\|\beta\|\) pour \(P e_{K} \geqslant 1\) et \(\delta\) suffisamment petit,

\[\left\|u-u_{h}\right\|_{G L S} \leqslant C h^{k+\frac{1}{2}}\|u\|_{H^{k+1}}\]

Convergence

En prenant \(\delta_{K}=\left(\|\beta\|_{\infty} / h_{K}+\varepsilon / h_{K}^{2}\right)^{-1},\) on obtient la convergence:

\[\left\|u-u_{h}\right\|_{L^{2}} \leqslant C h^{k+\frac{1}{2}}\|u\|_{H^{k+1}}\]

\(\rightarrow\) convergence sous-optimale d’un facteur \(1 / 2\)

Méthode 2: Méthode SUPG

SUPG: Streamline upwind Petrov-Galerkin
Definition: (Partie symétrique et antisymétrique)

\[\begin{array}{c} L u=-\varepsilon \Delta u+\nabla \cdot(\beta u)+\mu u \qquad=\underbrace{-\varepsilon \Delta u+\left[\mu+\frac{1}{2} \nabla \cdot \beta\right] u}_{=L_{S} u} \quad+\underbrace{\frac{1}{2}[\nabla \cdot(\beta u)+\beta \cdot \nabla u]}_{L_{S S} u} \\ L_{S}=\left(L+L^{T}\right) / 2: \text { symmetric part } \quad L_{S S}=\left(L-L^{T}\right) / 2: \text { skew-symmetric part } \end{array}\]

Remarques

\(\rightarrow \tilde{\mu}=\mu+\frac{1}{2} \nabla \cdot \beta\) intervient dans la partie symétrique de \(L\)
\(\rightarrow L_{S S} u=\beta \cdot \nabla u+\frac{1}{2} u \nabla \cdot \beta:\) si \(\beta\) est à divergence nulle, \(L_{S S}\) est réduit à la dérivée suivant \(\beta\).

Formulation variationnelle SUPG:

Trouver \(u \in V_{h},\)

\[a(u, v)+\underbrace{\sum_{K \in \mathcal{T}_{h}}\left(L u_{\mid K}, \delta_{K} L_{S S} v_{\mid K}\right)}_{a_{h}(u, v)}=\ell(v)+\underbrace{\sum_{K \in \mathcal{T}_{h}}\left(f, \delta_{K} L_{S S} v_{\mid K}\right)}_{\ell_{h}(v)}, \quad \forall v \in V_{h}\]

Remarques

\(\rightarrow\) peut s’obtenir à partir de
\[\left(L u, v+\delta L_{S S} v\right)=\left(f, v+\delta L_{S S} v\right)\]
\(\rightarrow L_{S S} v\) contient la dérivée suivant \(\beta \cdot \nabla v=\partial_{\beta} v\) d’où le nom streamline upwind
\(\rightarrow\) Petrov-Galerkin: les fonctions tests sont les \(v+\delta L_{S S} v\) avec \(v \in V_{h}\) au lieu de \(v \in V_{h}\)