Equation de diffusion

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$

et

$\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

Problème de diffusion non linéaire

Equation de Poisson non-linéaire

Verification

Formulation variationelle

Discrete Problem

Stratégies de résolution

Picard

Picard (2)

Picard (Iterations)

Newton algébrique 1

Newton algébrique 2

Newton algébrique 3

Newton algébrique 4

Newton algébrique 5

Newton EDP (1)

Newton EDP (2)

Newton EDP (3)

Newton EDP (4)

Equation de diffusion dépendant du temps

Equation de diffusion dépendant du temps

Formulation

Approche

Semi-discrétisation en temps

Résolution en espace

Remarques

Advection Diffusion Reaction

Problèmes type

Problèmes type

Advection diffusion stationnaire

Résolution

Problème

Nombre de Péclet

Méthode 1: Galerkin Least Square (GaLS)

Galerkin Least Square (GaLS)

Convergence

Convergence

Méthode 2: Méthode SUPG

Formulation variationnelle SUPG: