-- Simulation numérique en mécanique des fluides --
TP n°4
Dérivation numérique
simon.marie@lecnam.net
import numpy as np
import matplotlib.pyplot as plt
# la librairie misc de scipy permet de convertir les images en objets numpy
import scipy.misc as scm
# la librairie time nous sera utile pour calculer le temps d'execution des scripts
import time
#Option pour afficher les figures dans le notebook et eviter le plt.show():
%matplotlib inline
1 - Derivation numérique d'une fonction analytique
On se propose ici de tester un schéma aux différences finies pour calculer la dérivée d'une fonction sinus $f(x)=\sin(kx)$ soit (comme chacun sait) $f'(x)=k\cos(kx)$:
# On définie une première fonction qui va calculer une approximation numérique de la dérivée avec un schéma
# décentré d'ordre 1
def diff(f,dx):
df=np.zeros(f.shape)
# on peut l'écrire avec une boucle for
for i in range(len(f)-2):
df[i]= (f[i+1]-f[i])/dx # ordre 1 avec boucle
#df[0:-1]=(f[1:]-f[0:-1])/dx # ordre 1 avec slicing
#df[0:-2]=... # ordre 2
return df
# On définie une autre fonction qui va calculer une approximation numérique de la dérivée avec un schéma
# centré d'ordre 2
def diffc(f,dx):
df=np.zeros(f.shape)
for i in range(1,len(f)-1):
df[i]= (f[i+1]-f[i-1])/(2.*dx) # ordre 1
#df[1:-1]=... # ordre 2
return df
dx=0.05
x=np.arange(-1,1,dx)
y=np.sin(10.*x)
dy_th=10.*np.cos(10.*x)
dy=diff(y,dx)
dy2=diffc(y,dx)
plt.figure(figsize=(10,5))
plt.plot(x,y,'k',label='$f(x)$')
plt.plot(x,dy,'ro',label='ordre $1$ décentré')
plt.plot(x,dy2,'gs',label='ordre $2$ centré')
plt.plot(x,dy_th,'b',label='$f\'(x)$')
plt.legend()
Le schéma d'ordre 2 est plus précis lorsque $\Delta x$ est faible. On pourra essayer de changer les valeurs de $\Delta x$.
2 - Dérivation 2D
On se propose dans cette partie de définir des fonctions permetant de calculer le gradient et le laplacien d'une quantité $f$ donnée, dépendant de deux variables $x$ et $y$. On rapelle:
$$ \mathbf{grad}.x= \left( \begin{array}{c} \dfrac{\partial f}{\partial x} \\ \dfrac{\partial f}{\partial y} \\ \end{array} \right) $$
$$ \Delta f = \dfrac{\partial^2 f}{\partial x^2}+\dfrac{\partial^2 f}{\partial y^2} $$
On utilisera dans ces fonctions les différents schémas d'espace vues en cours.
# On construit la fonction grad avec 3 argument, 1 obligatoire: f
# et deux facultatifs o et dx qui sont l'ordre souhaité et la taille de maille.
def grad(f,o=1,dx=1.):
#On note dfx la composante suivant x du gradient:
dfx=np.zeros((f.shape[0],f.shape[1])) #on initialise dfx a 0
#On note dfy la composante suivant y du gradient:
dfy=np.zeros((f.shape[0],f.shape[1])) #on initialise dfy a 0
# Calcul de dfx et dfy en fonction de l'ordre:
if o==1: # ordre 1 decentre
dfx[0:-1,:]=(f[1:,:]-f[0:-1,:])/dx
dfy[:,0:-1]=(f[:,1:]-f[:,0:-1])/dx
elif o==2: # ordre 2 centre
dfx[1:-1,:]=(f[2:,:]-f[0:-2,:])/dx
dfy[:,1:-1]=(f[:,2:]-f[:,0:-2])/dx
else:
print('o < 2 ')
return dfx,dfy
def lap(f,dx=1.):
Lx=np.zeros((f.shape[0],f.shape[1])) #on initialise Lx a 0
Ly=np.zeros((f.shape[0],f.shape[1])) #on initialise Ly a 0
Lx[1:-1,:]=(f[0:-2,:]-2.*f[1:-1,:]+f[2:,:])/dx
Ly[:,1:-1]=(f[:,0:-2]-2.*f[:,1:-1]+f[:,2:])/dx
return Lx+Ly
Interprétation visuelle
Pour bien illustrer la notion de dérivation numérique, on peut utiliser une image en niveau de gris et appliquer les opérateur de dérivation sur cette image. On pourra utiliser l' image suivante (attention à bien l'enregistrer en local si vous executez le notebook):
 |
 |
f=scm.imread('CNAM.jpg')
plt.imshow(f,cmap='gray')
print('min/max:',f.min(),f.max())
print('taille:',f.shape)
fig=plt.figure(figsize=(15,10))
dfx,dfy=grad(f,2)
plt.imshow(dfx+dfy,cmap='Greys');plt.colorbar()
fig=plt.figure(figsize=(10,5))
fig.add_subplot(121)
dfx,dfy=grad(f,2)
plt.imshow((dfx+dfy)[200:550,300:500],cmap='Greys',vmin=-100,vmax=100);plt.colorbar()
fig.add_subplot(122)
Lap=lap(f)
plt.imshow(Lap[200:550,300:500],cmap='Greys',vmin=-100,vmax=100);plt.colorbar()
def filt_lap(f,dx=1.):
filt=np.zeros((f.shape[0],f.shape[1])) #on initialise Lx a 0
filt[1:-1,1:-1]=0.5*(f[0:-2,0:-2]+f[2:,2:])/dx
return (filt)/2.
Interprétation en mécanique des fluides
En represnant le champ de pression et de vitesse du TP n°3, calculer les opérateurs différentiels suivant associés au champ de vitesse:
Le gradient de densité: $\mathbf{grad}(\rho)=\left( \begin{array}{c} \dfrac{\partial \rho}{\partial x} \\ \dfrac{\partial \rho}{\partial y}\end{array} \right.$
Le laplacien de densité: $\Delta \rho = \dfrac{\partial^2 \rho}{\partial x^2}+\dfrac{\partial^2 \rho}{\partial y^2}$
- La divergence de la vitesse: div($\mathbf{U})=\dfrac{\partial U_x}{\partial x} +\dfrac{\partial U_y}{\partial y} $
- La vorticité: $\omega_z=\begin{array}{c} \dfrac{\partial U_x}{\partial y}-\dfrac{\partial U_y}{\partial x} \end{array} $
nx,ny=600,600
rho=np.fromfile('CM_Rho.dat').reshape(nx,ny)
Ux=np.fromfile('CM_UX.dat').reshape(nx,ny)
Uy=np.fromfile('CM_UY.dat').reshape(nx,ny)
drhodx,drhody=grad(rho)
Lap_rho=lap(rho)
duxdx,duxdy=grad(Ux)
duydx,duydy=grad(Uy)
grad_rho=np.zeros((nx,ny,2))
grad_rho[:,:,0]=drhodx
grad_rho[:,:,1]=drhody
divU=duxdx+duydy
Vortz=duxdy-duydx
plt.pcolormesh(Vortz.T,shading='gouraud')
plt.xlim([250,350]);plt.ylim([250,350])
plt.colorbar()
Conclusion¶
La dérivation numérique permet d'estimer les dérivée d'une grandeur à l'aide d'approximation dont on peut maitriser et connaitre l'erreur. La dérivée en espace ($\partial / \partial x$) met en évidence les variations d'une grandeur dans l'espace. La combinaison de ces dérivée permet de représenter des grandeurs plus complexe comme la divergence ou le rotationel traduisant un comportement particulier (rotationel, incompressible) de l'écoulement.
from IPython.core.display import HTML
style=open('../notebooks.css', "r").read()
HTML(style)