-- Mécanique des fluides numérique --

TP n°7

Résolution numérique de l'équation de Poisson

simon.marie@lecnam.net
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import copy
import time
%matplotlib inline

On se propose dans les 2 TP suivant d'étudier la résolution numérique de l'équation de Poisson.

Le TP7 est dédiée à la comparaison des méthodes directes et itératives pour trouver un état d'équilibre.

Le TP8 se concentrera sur l'utilisation de l'équation de Poisson pour calculer un champ de pression à partir d'un champ de vitesse donné.

Comparaison des méthodes de résolution d'un système elliptique

L'équation de Poisson s'écrit en 1D:

$$ \displaystyle{\dfrac{\partial^2 p}{\partial x^2}=f(x)} $$

à laquelle on associe pour ce TP les conditions aux limites suivantes:

$$ p(x=0)=2\\ \left.\dfrac{\partial p}{\partial x}\right|_{x=L}=0 $$

On choisira pour cette première partie: $$ f(x)=0.2(1+x^2-2x^3) $$

Résoudre ce problème revient à chercher une solution d'équilibre entre les contraintes imposées par le terme source $f(x)$ et les conditions aux limites.

Méthodes itératives

Une première façon de résoudre ce type de problème consiste à écrire l'équation discrète que l'on peut écrire en utilisant un schéma centré d'ordre 2 en espace:

$$ \dfrac{p_{i+1}-2p_{i}+p_{i-1}}{\Delta x^2}=f_{i} $$

ce qui donne pour $p_{i}$:

$$ \displaystyle{p_{i}=\dfrac{1}{2}(p_{i+1}+p_{i-1}-\Delta x^2 f_{i})} $$

Pour résoudre numériquement cette équation on défini $f_{i}$ comme étant une valeur d'entrée et on met à jour la valeur de $p_{i}$ jusqu'à satisfaire un critère de convergence de type $max|p_{i}^{nouveau}-p_{i}^{ancien}|\leq q$.

Définir une fonction permettant la résolution itérative de l'équation de Poisson ayant pour paramètres $q$, $\Delta x$ et $f_{i}$

On propose le canevas suivant:

In [ ]:
def poisson_it(f,q,dx,nmax=1000000):
    n=0
    p=np.ones(len(f))
    ... # CL type Dirichlet
    norm=q+1.
    while (norm>q):
        n+=1;pn=p.copy()
        p[1:-1]=...
        p[-1] = ... # CL type Neumann
        norm=...
        if (n>nmax):
            print('nmax atteint !!!!  q=',norm)
            break
    print(n)
    return(p)

Utiliser cette fonction pour résoudre l'équation de Poisson sur le segment $[-1,1]$. On prendra $\Delta x=0.01$

In [ ]:
dx=...
x=...
f=0.2*(1.+x**2-2.*x**3)

p=poisson_it(f,...,dx)
plt.plot(p)

Faire varier la valeur de $q$ et commenter. (On pourra prendre par exemple $q=10^{-5},q=10^{-6},q=10^{-7}$)

Discussions:

Méthodes directes

Un autre méthode pour résoudre les systèmes elliptiques consiste a résoudre le système linéaire en inversant direcement la matrice du système. En effet, on peut écrire le problème sous la forme matricielle: $$ \mathbf{A} \mathbf{p}=\mathbf{F} $$

avec:

$$ \mathbf{p}=(p_1,p_2,...,p_{nx})^T $$$$ \mathbf{F}=(f_1,f_2,...,f_{nx})^T $$

et: $$ A= \left[ \begin{array}{ccccc} -2&1&0&\cdots&0\\ 1&-2&1&\ddots&\vdots\\ 0&1&\ddots&1&0\\ \vdots&\ddots&1&-2&1\\ 0&\cdots&0&1&-2 \end{array} \right] $$

Définir une fonction permettant la résolution directe de l'équation de Poisson ayant pour paramètres $\Delta x$ et $f_{i}$. On utilisera les fonctions eye et diag de numpy pour la construction de la matrice ainsi que la fonction linalg.solve de numpy pour la résolution du système.

Pour les conditions aux limites:

On pourra écrire l'algorithme itératif aux deux extrémités:

  • En x=1 on a $p_{0}$ imposé: $$ \dfrac{p_{2}-2p_{1}+p_{0}}{\Delta x^2}=f_{1} $$

  • En x=L-1 on a $p_{L}=p_{L-1}$: $$ \dfrac{p_{L}-2p_{L-1}+p_{L-2}}{\Delta x^2}=\dfrac{-p_{L-1}+p_{L-2}}{\Delta x^2}=f_{L-1} $$

Proposer une solution pour incorporer les conditions aux limites dans la matrice $A$

In [ ]:
def poisson_dir(f,dx):
    A=...
    ...
    ... # CL
    ...
    p=np.linalg.solve(A,...)
    return p

Résoudre numériquement l'équation de Poisson sur le même domaine et comparer les solutions obtenues avec les deux méthodes. (On pourra comparer les temps d'execution avec la fonction time)

In [ ]:
pn=poisson_dir(f,dx)
plt.plot(x,pn,'k')
plt.plot(x[::10],p[::10],'or')

Discussions:

Conclusion

Présentez une synthèse de vos résultats et des étapes de calcul de votre TP.

Note: Vous pouvez télécharger ce notebook (fichier ipynb) et l'executer sur votre machine (installer Anaconda) ou sur le serveur JupyterHub du CNAM !

In [ ]:
from IPython.core.display import HTML
style=open('notebooks.css', "r").read()
HTML(style)