-- Simulation numérique en mécanique des fluides --
TP n°8
Résolution numérique de l'équation de Poisson 2D
éléments de correction
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 imcompressibles
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)=U_0(y-y_0)\exp{\left(-\dfrac{\ln{2}}{L^2}\left[(x-x_0)^2+(y-y_0)^2\right]\right)} \\ u_y(x,y)=-U_0(x-x_0)\exp{\left(-\dfrac{\ln{2}}{L^2}\left[(x-x_0)^2+(y-y_0)^2\right]\right)} $$On peut définir une fonction qui représente le champ de vitesse du tourbillon:
def tourbillon(nx,ny,L,U0):
x,y=np.arange(nx),np.arange(ny)
X=np.exp((-0.6931471806/L**2)*((x-float(nx/2.))**2))
Y=np.exp((-0.6931471806/L**2)*((y-float(ny/2.))**2))
X1=X*(x-float(nx/2.))
Y1=Y*(y-float(ny/2.))
ux=U0*np.einsum('i,j->ij',X,Y1)
uy=-U0*np.einsum('i,j->ij',X1,Y)
return ux,uy
nx,ny=100,100
p=np.zeros((nx,ny))
ux,uy=tourbillon(nx,ny,10,0.1)
x,y=np.arange(nx),np.arange(ny)
X,Y=np.meshgrid(x,y)
fig=plt.figure(figsize=(15,5))
fig.add_subplot(121)
plt.contour((ux**2+uy**2)**0.5);plt.axis('equal')
plt.colorbar()
fig.add_subplot(122)
plt.streamplot(X,Y,ux.T,uy.T);plt.axis('equal')
#plt.xlim([20,80]);plt.ylim([20,80])
On observe un champ de vitesse caractéristique d'un tourbillon stationnaire, avec une vitesse nulle au centre et une vitesse max annulaire.
Calculon à présent la source engendré par ce type de champ:
def source(ux,uy,dx):
nx,ny=ux.shape
f=np.zeros((nx,ny))
f[1:-1,1:-1]=0.25*(ux[2:,1:-1]-ux[0:-2,1:-1])**2/dx**2 # (dux/dx)^2
f[1:-1,1:-1]+=0.5*(ux[1:-1,2:]-ux[1:-1,0:-2])*(uy[2:,1:-1]-uy[0:-2,1:-1])/dx**2 # 2(dux/dy)(duy/dx)
f[1:-1,1:-1]+=0.25*(uy[1:-1,2:]-uy[1:-1,0:-2])**2/dx**2 # (duy/dy)^2
return f
f=-source(ux,uy,0.1)
plt.pcolormesh(f,cmap="jet",shading='gouraud');plt.axis('equal')
plt.colorbar()
On observe une source non nulle à l'endroit du tourbillon, pariculièrement important là ou les gradient de vitesses sont important. Le champ de pression représenté ici ne correspond pas à celui associé à la présence d'un tourbillon (dépression).
On va donc résoudre l'équation de Poisson en 2 dimensions, en utilisant une méthode itérative:
Calcul du champ de pression associé
Résolution itérative
L'algorithme itératif pour ce type de problème s'écrit (pour un schéma centré d'ordre 2):
$$ \dfrac{p_{i+1,j}-2p_{i,j}+p_{i-1,j}}{dx^2}+\dfrac{p_{i,j+1}-2p_{i,j}+p_{i,j-1}}{dy^2}=f_{ij} $$soit (en prenant $dx=dy$):
$$ p{i,j}=\dfrac{1}{4}(p_{i+1,j}+p_{i-1,j}+p_{i,j+1}+p_{i,j-1}-dx^2 f_{i,j}) $$Pour résoudre ce problème on va choisir ici d'imposer des gradients de pression nulle sur tous les bords du domaine. En effet, le but est d'obtenir le champ de pression associé à la présence d'un tourbillon. Ici les bords du domaine sont très peu impacté par le tourbillon. Le choix des conditions au bords importe donc peu. (On pourra aussi utiliser Dirichlet au bord en imposant une pression uniforme $p=1$, ça ne changera rien).
def poisson_2D(p,f,q,dx,nmax=10000):
n=0;conv=1;norm=1.
eps=0.1
while (norm > q):
n+=1;pn=p.copy()
p[1:-1,1:-1]=0.25*(pn[0:-2,1:-1]+pn[2:,1:-1]+pn[1:-1,2:]+pn[1:-1,0:-2]-f[1:-1,1:-1]*dx**2)
p[:,-1] = p[:,-2] # dp/dy = 0 at y=end
p[:,0] = p[:,1] # dp/dy = 0 at y = 0
p[0,:] = p[1,:] # dp/dx = 0 at x = 0
p[-1,:] = p[-2,:] # dp/dx = 0 at x = end
norm=abs(p-pn).max()
if (n>nmax):
print('nmax atteint: q='+str(norm))
conv=0
break
print(n)
Résolvons maintenant Poisson pour un champ de vitesse tourbillonaire:
poisson_2D(p,f,1e-7,0.1)
fig=plt.figure(figsize=(12,5),dpi=120)
fig.add_subplot(121)
plt.pcolormesh(uy,cmap="jet",shading='gouraud');plt.axis('equal')
plt.colorbar()
fig.add_subplot(122)
plt.pcolormesh(p,cmap="jet",shading='gouraud');plt.axis('equal')
plt.colorbar()
On obtient donc un champ de pression cohérent avec la présence d'un tourbillon. La pression baisse au centre du tourbillon qui est donc bien caractérisé par une depression:
plt.plot(p[50,:])
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)^2$ 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.
Commencons par calculer la matrice de façon intuitive en remplissant les diagonales comme dans le TP7:
def Lap(nx,ny):
A=-4.*np.eye(nx*ny)
A+=np.diag(np.ones(nx*ny-1),1)
A+=np.diag(np.ones(nx*ny-1),-1)
A+=np.diag(np.ones(nx*ny-ny),ny)
A+=np.diag(np.ones(nx*ny-ny),-ny)
return A
nx,ny,dx=100,100,0.1
ux,uy=tourbillon(nx,ny,10,0.1)
f=-source(ux,uy,dx).reshape(nx*ny)
On résout alors le système linéaire en mesurant le temps de calcul. On inclu dans le temps de calcul la création de la matrice A. Nous verrons que cela est important dans la suite.
t0=time.time()
A=Lap(nx,ny)
pd=np.linalg.solve(A,dx**2*f)
time.time()-t0
C'est relativement long ! En fait il faut créer une très grosse matrice en mémoire et ensuite résoudre le système linéaire associé ! Ce qui fait beaucoup de temps pour finalement une matrice presque vide !!!
plt.pcolormesh(pd.reshape(nx,ny),cmap='jet',shading='gouraud');plt.axis('equal')
plt.colorbar()
En fait, il existe des outils pour manipuler les matrices creuses. On utilisera ici les outils de la bibliothèque scipy.sparse:
from scipy.sparse import dia_matrix
from scipy.sparse.linalg import spsolve
L'idée de ces outils c'est de créer uniquement les diagonales en mémoire et d'utiliser des algorithmes adaptés aux matrices creuses. On pourra regarder par exemple la documentation de spsolve
def Lap_sp(nx,ny):
data = np.array([-4*np.ones(nx*ny),np.ones(nx*ny),np.ones(nx*ny),np.ones(nx*ny),np.ones(nx*ny)])
offsets = np.array([0, -1, 1,ny,-ny])
A=dia_matrix((data, offsets), shape=(nx*ny, nx*ny))
return A
t0=time.time()
B=Lap_sp(nx,ny)
pd_s = spsolve(B, dx**2*f)
time.time()-t0
C'est beaucoup plus rapide !!!! Nous reviendrons plus largement sur ce type de technique dans le prochain TP.
plt.pcolormesh(pd_s.reshape(nx,ny),cmap='jet',shading='gouraud');plt.axis('equal')
plt.colorbar()
fig=plt.figure(figsize=(12,5),dpi=120)
plt.plot(x[::2],p[50,::2],'s',label='Methode itérative')
plt.plot(pd.reshape(nx,ny)[50,:],label='Methode directe')
plt.plot(x[::4],pd_s.reshape(nx,ny)[50,::4],'o',label='Methode directe (sparse)')
plt.legend()
Conclusion¶
Ce TP doit vous avoir permis de mieux comprendre le rôle de l'équation de Poisson pour la pression pour les écoulements incompressibles. Elle permet de déduire le champ de pression à partir du champ de vitesse. Vous avez ainsi pu retrouver le champ de pression associé à un tourbillon stationnaire.
Ce TP a également permis d'appliquer les solveurs utilisé sur un cas 2D. En particulier, vous avez vu que la résolution directe du problème à partir de la matrice du Laplacien était couteurse en ressources. Cependant, étant donnée la forme particulière de cette matrice, il est plus judicieux d'utiliser les outils "sparse" adaptés aux matrices creuse dans Python.
Dans le prochain TP, il sera question de la Méthode de Newton que nous utiliserons pour chercher une solution stationnaire de Navier Stokes. De nombreux concepts introduit ici seron réutilisé.
from IPython.core.display import HTML
style=open('../notebooks.css', "r").read()
HTML(style)