-- Mécanique des fluides numérique --
TP n°8
Résolution numérique de l'équation de Poisson 2D
simon.marie@lecnam.net
import numpy as np
import matplotlib.pyplot as plt
import copy
import time
%matplotlib inline
L'équation de Poisson pour les écoulements incompressibles
Pour la simulation numérique des écoulements en régime incompressible, il faut trouver un moyen de relier la pression et la vitesse. C'est l'équation de Poisson qui joue ce rôle ! Pour expliciter le terme source $f$ on calcul la divergence de l'équation de quantité de mouvement: $\partial_x(Qdm_x)+\partial_y(Qdm_y)$ et en utilisant l'hypothèse d'incompressibilité $\partial_x u_x+\partial_y u_y=0$, on obtient l'équation de Poisson avec le terme source suivant:
$$ \displaystyle{-\dfrac{1}{\rho}\left(\dfrac{\partial^2 p}{\partial x^2}+\dfrac{\partial^2 p}{\partial y^2}\right)=\left(\dfrac{\partial u_x}{\partial x}\right)^2+2\dfrac{\partial u_x}{\partial y}\dfrac{\partial u_y}{\partial x}+\left(\dfrac{\partial u_y}{\partial y}\right)^2} $$
On va donc dans cette partie, utiliser l'équation précédente pour calculer un champ de pression à partir d'un champ de vitesse connue.
Tourbillon stationnaire
Soit le champ de vitesse $(u_x,u_y)$ défini par:
$$ u_x(x,y)=(y-y_0)U_0\exp{\left(-\dfrac{\ln{2}}{L^2}\left[(x-x_0)^2+(y-y_0)^2\right]\right)} \\ u_y(x,y)=-(x-x_0)U_0\exp{\left(-\dfrac{\ln{2}}{L^2}\left[(x-x_0)^2+(y-y_0)^2\right]\right)} $$
Définir une fonction permettant de calculer ce champ de vitesse à partir des vecteurs $x$ et $y$ et des paramètres $L$ et $U_0$. On prendra le point $(x_0,y_0)$ au centre du domaine.
def tourbillon(x,y,L,U0):
nx=len(x)
ny=len(y)
...
...
ux=...
uy=...
return ux,uy
Tracer le champ de vitesse dans un domaine $([-100;100],[-100;100])$ avec $L=10$, $U_0=0.1$ et $\Delta x=0.1$
dx=0.1
x=np.arange(-100,-100,dx)
y=np.arange(-100,-100,dx)
ux,uy=tourbillon(x,y,10,0.1)
plt.pcolormesh((ux**2+uy**2)**0.5,shading='gouraud');plt.axis('equal')
plt.colorbar()
On pourra également utiliser la fonction streamplot pour visualiser les lignes de courant:
plt.streamplot(X,Y,ux.T,uy.T);plt.axis('equal')
Définir une fonction permettant le calcul du terme source à partir du champ de vitesse $(u_x,u_y)$ et tracer le terme source
def source(ux,uy,dx):
nx,ny=ux.shape
f=np.zeros((nx,ny))
f[1:-1,1:-1]=...
return f
f=...
plt.pcolormesh(f,shading='gouraud');plt.axis('equal')
plt.colorbar()
Discussion:
Calcul du champ de pression associé
1 - Méthode itérative
Une résolution itérative du problème 2D peut s'écrire de façon analogue à la première partie:
$$ \dfrac{p_{i+1,j}-2p_{i,j}+p_{i-1,j}}{\Delta x^2}+\dfrac{p_{i,j+1}-2p_{i,j}+p_{i,j-1}}{\Delta y^2}=f_{i,j} $$
On choisit d'imposer comme condition aux limites un gradient nulle à chaque bord.
Sur la base de ce que vous avez fait dans le TP7, proposer une fonction pour la résolution de l'équation de Poisson en 2D
def poisson_2D(f,q,dx,nmax=10000):
n=0;norm=1.
p=...
while (norm > q):
n+=1;pn=p.copy()
p[1:-1,1:-1]=...
p[:,-1] = ...
p[:,0] = ...
p[0,:] = ...
p[-1,:] = ...
norm=abs(p-pn).max()
if (n>nmax):
print()'nmax atteint: q='+str(norm))
break
return p
Calculer alors le champ de pression associé au tourbillon stationnaire
...
...
...
On pourra faire varier les critères de convergence et discuter des résultats
Discussions:
2 - Résolution Matricielle
On souhaite maintenant, comme dans le TP7, résoudre le même problème mais avec une méthode directe. Ceci implique tout d'abord d'écrire le problème sous forme matricielle du type:
$$ \mathbf{A}P=F $$
Que l'on écrira:
$$ \mathbf{P}=(p_{11},p_{12},...,p_{1,ny},p_{21},...,...,p_{nx,ny})^T $$
$$ \mathbf{F}=(f_{11},f_{12},...,f_{1,ny},f_{21},...,...,f_{nx,ny})^T $$
et: $$ A= \left[ \begin{array}{ccccc} -4&1&0&\cdots&1&0&\cdots&0\\ 1&-4&1&\ddots&\vdots&\ddots&\ddots&\vdots\\ 0&1&\ddots&1&0&\vdots&\ddots&0\\ \vdots&\ddots&1&-4&\ddots&\ddots&0&1\\ 1&\ddots&0&\ddots&-4&1&0&\ddots\\ 0&\ddots&0&0&1&\ddots&1&0\\ \vdots&\ddots&\ddots&0&0&1&-4&1\\ 0&\cdots&0&1&0&0&1&-4\\ \end{array} \right] $$
Ainsi la matrice $A$ est une matrice de taille $n_x \times n_y$ de type "penta-diagonal" c'est à dire avec 5 diagonales non nulle. La matrice $A$ est donc relativement "vide" et contient beaucoup plus de valeur nulle que de valeur non nulles. On apelle cela une matrice "creuse" ou "sparse" en anglais.
nx,ny,dx=100,100,0.1
ux,uy=tourbillon(nx,ny,10,0.1)
f=-source(ux,uy,dx).reshape(nx*ny)
Surla base de ce que vous avez fait dans le TP7, définir une fonction permettant de remplir la matrice A associé au laplacien:
def Lap(nx,ny):
A=...
return A
Résoudre le système et mesurer le temps nécessaire.
t0=time.time()
A=Lap(nx,ny)
pd=np.linalg.solve(...,...)
...
Commentaires:....
En pratique on pourra regarder les outils permettant de manipuler les matrices creuses:
from scipy.sparse import dia_matrix
from scipy.sparse.linalg import spsolve
def Lap_sp(nx,ny):
data = ...
offsets = ...
A=dia_matrix((data, offsets), shape=(nx*ny, nx*ny))
return A
t0=time.time()
B=Lap_sp(nx,ny)
pd_s = spsolve(..., ...)
time.time()-t0
Comparer les 3 résultats:
fig=plt.figure(figsize=(12,5),dpi=120)
plt.plot(...,label='Methode itérative')
plt.plot(...,label='Methode directe')
plt.plot(...,label='Methode directe (sparse)')
plt.legend()
Commentaires.
Conclusion¶
Présentez une synthèse de vos résultats et des étapes de calcul de votre TP.
from IPython.core.display import HTML
style=open('notebooks.css', "r").read()
HTML(style)