Aérodynamisme : étude de profils et comparaisons

                         Computational Fluid Dynamics (CFD) applied to Formula 1 study
                     "Dynamique des fluides numérique appliquée à l'étude d'une Formule 1"

Informations générales

1)
(pn[1:-1, 2:] + pn[1:-1, 0:-2]) * dy2 + (pn[2:, 1:-1] + pn[0:-2, 1:-1]) * dx2) /
                      (2 * (dx**2 + dy**2)) -
                      dx**2 * dy**2 / (2 * (dx**2 + dy**2)) * (rho * 
                      (1 / dt * ((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx) +
                      (v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy)) -
                      ((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx))**2 -
                      2 * ((u[2:, 1:-1] - u[0:-2, 1:-1]) / (2 * dy) *
                      (v[1:-1, 2:] - v[1:-1, 0:-2]) / (2 * dx))-
                      ((v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy))**2)))
                          
      #conditions aux limites en xmin
      p[1:-1, 0] = (((pn[1:-1, 1] + pn[1:-1, -1])* dy**2 +
                   (pn[2:, 0] + pn[0:-2, 0]) * dx**2) /
                   (2 * (dx**2 + dy**2)) -
                   dx**2 * dy**2 / (2 * (dx**2 + dy**2)) * 
                   (rho * (1 / dt * ((u[1:-1, 1] - u[1:-1, -1]) / (2 * dx) +
                   (v[2:, 0] - v[0:-2, 0]) / (2 * dy)) -
                   ((u[1:-1, 1] - u[1:-1, -1]) / (2 * dx))**2 -
                   2 * ((u[2:, 0] - u[0:-2, 0]) / (2 * dy) *
                   (v[1:-1, 1] - v[1:-1, -1]) / (2 * dx))-
                   ((v[2:, 0] - v[0:-2, 0]) / (2 * dy))**2)))
      #conditions aux limites en xmax
      p[1:-1, -1] = (((pn[1:-1, 0] + pn[1:-1, -2])* dy**2 +
                    (pn[2:, -1] + pn[0:-2, -1]) * dx**2) /
                    (2 * (dx**2 + dy**2)) -
                    dx**2 * dy**2 / (2 * (dx**2 + dy**2)) * (rho * (1 / dt *
                    ((u[1:-1, 0] - u[1:-1,-2]) / (2 * dx) +
                    (v[2:, -1] - v[0:-2, -1]) / (2 * dy)) -
                    ((u[1:-1, 0] - u[1:-1, -2]) / (2 * dx))**2 -
                    2 * ((u[2:, -1] - u[0:-2, -1]) / (2 * dy) *
                    (v[1:-1, 0] - v[1:-1, -2]) / (2 * dx)) -
                    ((v[2:, -1] - v[0:-2, -1]) / (2 * dy))**2)))
      #conditions aux limites sur y
      p[0,:]=p[1,:]  #dérivée selon y nulle à ymin
      p[-1,:]=p[-2,:]  #dérivée selon y nulle à ymax
  
  return p
</code> Le travail numérique revient ainsi à préciser la pression et les vitesses en tout point d'un espace discret. Cependant, en plus du travail sur les équations de Navier-Stokes, nous allons devoir préciser les forces aérodynamiques en tout point (dont les vitesses dépendent) ainsi que sur les surfaces de références $S_{ref}$ dont dépendent les valeurs des forces. On cherche à visualiser les effets des flux d'air sur un cube (un carré en deux dimensions) et sur une sphère (un disque en deux dimensions): naturellement, on s'attend à ce que les courants évitent les objets et les contournent. Définissons maintenant la force aérodynamique sur le modèle de notre formule par boucles successives. 
#définition de la force aérodynamique
def force(u,v,rho,s,cx,cy):
    fx=np.empty(np.shape(u))
    fy=np.empty(np.shape(v))
    lin_x,col_x=np.shape(fx)
    lin_y,col_y=np.shape(fy)

    for i in range(lin_x):
        for j in range(col_x):
            fx[i,j]=0.5*rho*s[i,j]*cx*(u[i,j])**2
    for i in range(lin_y):
        for j in range(col_y):
            fy[i,j]=0.5*rho*s[i,j]*cy*(v[i,j])**2
              
    return (fx,fy)
Le traçage de notre surface de référence est un peu plus compliquée. De manière évidente, nous ne pourrons intégrer des surfaces complexes à notre modèle tel une Formule 1, des voitures ou encore des ailes profilées. La définition de nos surfaces de références revient à considérer une matrice de surface de taille $n_x \times n_y$: à chaque valeur de notre matrice correspond une valeur booléenne ($0$ si aucune surface n'est présente localement, $1$ dans le cas contraire). Ainsi, nous nous servons de notre grille pour imprimer notre surface, ce qui explique que les profils incurvés comme le cercle sont définis en forme de “dents” ou encore en “escaliers”. Cependant, pour contrer cet aspect cassant, il suffit d'augmenter la densité de nos points c'est à dire d'augmenter la taille de nos matrices. Le problème majeur est que cela revient à calculer des matrices de plus en plus grandes: nous devons donc faire un choix entre la résolution de notre figure final et le temps, les ressources de l’exécution. Pour les figures $(30,30)$, le programme finit en moins d'une minute mais les surfaces sont rugueuses et nous avons peu de vecteurs dans notre champ. Pour information, le résultat avec les matrices $(100,100)$ donne plus de précision mais se termine en cinq bonnes minutes ! 
#définition de la surface de référence
def surface(x,y,cmX,cmY,taille,chaine):
    s=np.empty((np.size(y),np.size(x)))
    if chaine=='cube':
        xm=cmX-(taille/2)
        xM=cmX+(taille/2)
        ym=cmY-(taille/2)
        yM=cmY+(taille/2)
        for i in range(np.size(x)):
            for j in range(np.size(y)):
                if x[i]<=xM and x[i]>=xm and y[j]<=yM and y[j]>=ym:
                    s[i,j]=1
                else:
                    s[i,j]=0
    if chaine=='cercle':
        for i in range(np.size(x)):
            for j in range(np.size(y)):
                if ((((x[i]-cmX)**2)+((y[j]-cmY)**2))<=taille**2):
                    s[i,j]=1
                else:
                    s[i,j]=0
            
    return s
Initialement, on suppose un champ de vitesse selon $\vec{e_x}$ uniquement (dans le cas de la soufflerie, cette hypothèse est vérifiée, l'air est éjecté massivement vers l'avant selon la direction primaire et peut sur les côtés). Nous allons donc devoir réaliser une boucle tournant sur un critère de précision (plus précisément l'écart successif entre deux matrices de vitesses) et placer un compteur. 
#implementation
#variable de contrôle
accu=0.00001
#index de boucle
diff_u=1.0
stepcount=0

while diff_u>=accu:
    un=u.copy() #m.s-1
    vn=v.copy() #m.s-1
    p=pressure(p,u,v,rho,dx,dy,dt,n_time) #Pa
    fx,fy=force(u,v,rho,s,cx,cy) #N

    #intérieur des tableaux
    u[1:-1, 1:-1] = (un[1:-1, 1:-1] -
                    un[1:-1, 1:-1] * dt / dx * 
                    (un[1:-1, 1:-1] - un[1:-1, 0:-2]) -
                    vn[1:-1, 1:-1] * dt / dy * 
                    (un[1:-1, 1:-1] - un[0:-2, 1:-1]) -
                    dt / (2 * rho * dx) * 
                    (p[1:-1, 2:] - p[1:-1, 0:-2]) +
                    nu * (dt / dx**2 * 
                    (un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) +
                    dt / dy**2 * 
                    (un[2:, 1:-1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1])) +
                    fx[1:-1, 1:-1] * dt)

    v[1:-1, 1:-1] = (vn[1:-1, 1:-1] -
                     un[1:-1, 1:-1] * dt / dx * 
                    (vn[1:-1, 1:-1] - vn[1:-1, 0:-2]) -
                     vn[1:-1, 1:-1] * dt / dy * 
                    (vn[1:-1, 1:-1] - vn[0:-2, 1:-1]) -
                     dt / (2 * rho * dy) * 
                    (p[2:, 1:-1] - p[0:-2, 1:-1]) +
                     nu * (dt / dx**2 *
                    (vn[1:-1, 2:] - 2 * vn[1:-1, 1:-1] + vn[1:-1, 0:-2]) +
                     dt / dy**2 * 
                    (vn[2:, 1:-1] - 2 * vn[1:-1, 1:-1] + vn[0:-2, 1:-1])) + 
                    fy[1:-1, 1:-1] * dt)

    #conditions aux limites en xmin
    u[1:-1, 0] = (un[1:-1, 0] - un[1:-1, 0] * dt / dx *
                 (un[1:-1, 0] - un[1:-1, -1]) -
                 vn[1:-1, 0] * dt / dy * 
                 (un[1:-1, 0] - un[0:-2, 0]) - 
                 dt / (2 * rho * dx) * 
                 (p[1:-1, 1] - p[1:-1, -1]) + 
                 nu * (dt / dx**2 * 
                 (un[1:-1, 1] - 2 * un[1:-1, 0] + un[1:-1, -1]) +
                 dt / dy**2 *
                 (un[2:, 0] - 2 * un[1:-1, 0] + un[0:-2, 0])) + 
                 fx[1:-1, 0] * dt)
                
    v[1:-1, 0] = (vn[1:-1, 0] - un[1:-1, 0] * dt / dx *
                 (vn[1:-1, 0] - vn[1:-1, -1]) -
                 vn[1:-1, 0] * dt / dy *
                 (vn[1:-1, 0] - vn[0:-2, 0]) -
                 dt / (2 * rho * dy) * 
                 (p[2:, 0] - p[0:-2, 0]) +
                 nu * (dt / dx**2 * 
                 (vn[1:-1, 1] - 2 * vn[1:-1, 0] + vn[1:-1, -1]) +
                 dt / dy**2 * 
                 (vn[2:, 0] - 2 * vn[1:-1, 0] + vn[0:-2, 0])) + 
                 fy[1:-1, 0] * dt)

    #conditions aux limites en xmax
    u[1:-1, -1] = (un[1:-1, -1] - un[1:-1, -1] * dt / dx * 
                  (un[1:-1, -1] - un[1:-1, -2]) -
                  vn[1:-1, -1] * dt / dy * 
                  (un[1:-1, -1] - un[0:-2, -1]) -
                  dt / (2 * rho * dx) *
                  (p[1:-1, 0] - p[1:-1, -2]) + 
                  nu * (dt / dx**2 * 
                  (un[1:-1, 0] - 2 * un[1:-1,-1] + un[1:-1, -2]) +
                  dt / dy**2 * 
                  (un[2:, -1] - 2 * un[1:-1, -1] + un[0:-2, -1])) + 
                  fx[1:-1, -1] * dt)
                   
    v[1:-1, -1] = (vn[1:-1, -1] - un[1:-1, -1] * dt / dx *
                  (vn[1:-1, -1] - vn[1:-1, -2]) - 
                  vn[1:-1, -1] * dt / dy *
                  (vn[1:-1, -1] - vn[0:-2, -1]) -
                  dt / (2 * rho * dy) * 
                  (p[2:, -1] - p[0:-2, -1]) +
                  nu * (dt / dx**2 *
                  (vn[1:-1, 0] - 2 * vn[1:-1, -1] + vn[1:-1, -2]) +
                  dt / dy**2 *
                  (vn[2:, -1] - 2 * vn[1:-1, -1] + vn[0:-2, -1])) + 
                  fy[1:-1, -1] * dt)
    
    #on élimine les mauvaises valeurs à l'intérieur des surfaces
    for i in range(nx):
        for j in range(ny):
            if s[i,j]!=0:
                u[i,j]=0
                v[i,j]=0
    
    #vitesses initiales uniquement selon x
    u[:,0]=u0
    v[:,0]=0
    
    #test et actualisation de la boucle
    diff_u=(np.sum(u)-np.sum(un))/np.sum(u)
    stepcount=stepcount+1
On obtient ainsi une représentation du champ de vecteurs en tout point de l'espace que l'on a discrétisé (avec toutefois un nombre important de valeurs pour avoir un bel effet de densité de points). On met ça en forme sur un graphique avec échelle de couleurs. 
#graphe
X,Y=np.meshgrid(x,y)
color=np.hypot(u,v) #méthode pour colorier les vecteurs
plt.quiver(X,Y,u,v,color)
plt.xlabel('X (m)')
plt.ylabel('Y (m)')
plt.title('['+stype.upper()+']'+'   Traînée: '+
          str(cx) + ' Dérive: '+str(cy)+' Compteur: '
          +str(stepcount),fontsize=10)
plt.suptitle('Représentation graphique des courants',fontweight='bold')
plt.colorbar(label='V (m/s)')
plt.scatter(X,Y,s,c='black',marker='s',alpha=0.5) #surface
plt.legend()
plt.savefig(f_name+'_'+stype+'_'+str(nx)+'.jpg',dpi=300)
plt.show()
Détail complet du code sans modifications 
# -*- coding: utf-8 -*-

#importation des bibliothèques
import numpy as np
import matplotlib.pyplot as plt

#nettoyage des anciennes fenêtres pyplot
plt.close('all')

#entrées de l'utilisateur
f_name='aerodynamics_cfd' #nom du fichier
stype='cube' #type de surface considéré entre cube et cercle
cx=0.47 #coefficient de traînée
cy=0.47 #coefficient de dérive
u0=0.1 #m.s-1, vitesse initiale
    
#définition de la pression de Poisson
def pressure(p,u,v,rho,dx,dy,dt,n_time):
    pn=np.empty(np.shape(p)) #on crée une matrice similaire à p
    for k in range(n_time):
        pn = p.copy() #copie p pour faire des opérations dessus
        
        #intérieur du tableau
        p[1:-1, 1:-1] = (((pn[1:-1, 2:] + pn[1:-1, 0:-2]) * dy**2 +
                        (pn[2:, 1:-1] + pn[0:-2, 1:-1]) * dx**2) /
                        (2 * (dx**2 + dy**2)) -
                        dx**2 * dy**2 / (2 * (dx**2 + dy**2)) * (rho * 
                        (1 / dt * ((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx) +
                        (v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy)) -
                        ((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx))**2 -
                        2 * ((u[2:, 1:-1] - u[0:-2, 1:-1]) / (2 * dy) *
                        (v[1:-1, 2:] - v[1:-1, 0:-2]) / (2 * dx))-
                        ((v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy))**2)))
                            
        #conditions aux limites en xmin
        p[1:-1, 0] = (((pn[1:-1, 1] + pn[1:-1, -1])* dy**2 +
                     (pn[2:, 0] + pn[0:-2, 0]) * dx**2) /
                     (2 * (dx**2 + dy**2)) -
                     dx**2 * dy**2 / (2 * (dx**2 + dy**2)) * 
                     (rho * (1 / dt * ((u[1:-1, 1] - u[1:-1, -1]) / (2 * dx) +
                     (v[2:, 0] - v[0:-2, 0]) / (2 * dy)) -
                     ((u[1:-1, 1] - u[1:-1, -1]) / (2 * dx))**2 -
                     2 * ((u[2:, 0] - u[0:-2, 0]) / (2 * dy) *
                     (v[1:-1, 1] - v[1:-1, -1]) / (2 * dx))-
                     ((v[2:, 0] - v[0:-2, 0]) / (2 * dy))**2)))

        #conditions aux limites en xmax
        p[1:-1, -1] = (((pn[1:-1, 0] + pn[1:-1, -2])* dy**2 +
                      (pn[2:, -1] + pn[0:-2, -1]) * dx**2) /
                      (2 * (dx**2 + dy**2)) -
                      dx**2 * dy**2 / (2 * (dx**2 + dy**2)) * (rho * (1 / dt *
                      ((u[1:-1, 0] - u[1:-1,-2]) / (2 * dx) +
                      (v[2:, -1] - v[0:-2, -1]) / (2 * dy)) -
                      ((u[1:-1, 0] - u[1:-1, -2]) / (2 * dx))**2 -
                      2 * ((u[2:, -1] - u[0:-2, -1]) / (2 * dy) *
                      (v[1:-1, 0] - v[1:-1, -2]) / (2 * dx)) -
                      ((v[2:, -1] - v[0:-2, -1]) / (2 * dy))**2)))

        #conditions aux limites sur y
        p[0,:]=p[1,:]  #dérivée selon y nulle à ymin
        p[-1,:]=p[-2,:]  #dérivée selon y nulle à ymax
    
    return p

#définition de la force aérodynamique
def force(u,v,rho,s,cx,cy):
    fx=np.empty(np.shape(u))
    fy=np.empty(np.shape(v))
    lin_x,col_x=np.shape(fx)
    lin_y,col_y=np.shape(fy)

    for i in range(lin_x):
        for j in range(col_x):
            fx[i,j]=0.5*rho*s[i,j]*cx*(u[i,j])**2
    for i in range(lin_y):
        for j in range(col_y):
            fy[i,j]=0.5*rho*s[i,j]*cy*(v[i,j])**2
              
    return (fx,fy)

#définition de la surface de référence
def surface(x,y,cmX,cmY,taille,chaine):
    s=np.empty((np.size(y),np.size(x)))
    if chaine=='cube':
        xm=cmX-(taille/2)
        xM=cmX+(taille/2)
        ym=cmY-(taille/2)
        yM=cmY+(taille/2)
        for i in range(np.size(x)):
            for j in range(np.size(y)):
                if x[i]<=xM and x[i]>=xm and y[j]<=yM and y[j]>=ym:
                    s[i,j]=1
                else:
                    s[i,j]=0
    if chaine=='cercle':
        for i in range(np.size(x)):
            for j in range(np.size(y)):
                if ((((x[i]-cmX)**2)+((y[j]-cmY)**2))<=taille**2):
                    s[i,j]=1
                else:
                    s[i,j]=0
            
    return s
            
#variables physiques
l=0.3 #m
rho=1.204 #kg.m-3, masse volumique
eta=0.1 #Pa.s, viscosité
nu=eta/rho #m2.s-1, viscosité cinématique

#paramètres d'intégration
n_time=30
nx,ny=30,30
xmin,xmax=0,1 #m
ymin,ymax=0,1 #m
x=np.linspace(xmin,xmax,nx) #m
y=np.linspace(ymin,ymax,ny) #m
dt=0.01
dx=0.1
dy=0.1

#surface
cmX=0.5 #m
cmY=0.5 #m
l=0.3 #m
s=surface(x,y,cmX,cmY,l,stype) #m2

#initialisation des matrices
u=np.zeros((ny,nx))
un=np.zeros((ny,nx))
v=np.zeros((ny,nx))
vn=np.zeros((ny,nx))
p=np.ones((ny,nx))
pn=np.ones((ny,nx))

#on informe l'utilisateur que le programme est lancé
print('process launched, please wait')

#implementation
#variable de contrôle
accu=0.00001
#index de boucle
diff_u=1.0
stepcount=0

while diff_u>=accu:
    un=u.copy() #m.s-1
    vn=v.copy() #m.s-1
    p=pressure(p,u,v,rho,dx,dy,dt,n_time) #Pa
    fx,fy=force(u,v,rho,s,cx,cy) #N

    #intérieur des tableaux
    u[1:-1, 1:-1] = (un[1:-1, 1:-1] -
                    un[1:-1, 1:-1] * dt / dx * 
                    (un[1:-1, 1:-1] - un[1:-1, 0:-2]) -
                    vn[1:-1, 1:-1] * dt / dy * 
                    (un[1:-1, 1:-1] - un[0:-2, 1:-1]) -
                    dt / (2 * rho * dx) * 
                    (p[1:-1, 2:] - p[1:-1, 0:-2]) +
                    nu * (dt / dx**2 * 
                    (un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) +
                    dt / dy**2 * 
                    (un[2:, 1:-1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1])) +
                    fx[1:-1, 1:-1] * dt)

    v[1:-1, 1:-1] = (vn[1:-1, 1:-1] -
                     un[1:-1, 1:-1] * dt / dx * 
                    (vn[1:-1, 1:-1] - vn[1:-1, 0:-2]) -
                     vn[1:-1, 1:-1] * dt / dy * 
                    (vn[1:-1, 1:-1] - vn[0:-2, 1:-1]) -
                     dt / (2 * rho * dy) * 
                    (p[2:, 1:-1] - p[0:-2, 1:-1]) +
                     nu * (dt / dx**2 *
                    (vn[1:-1, 2:] - 2 * vn[1:-1, 1:-1] + vn[1:-1, 0:-2]) +
                     dt / dy**2 * 
                    (vn[2:, 1:-1] - 2 * vn[1:-1, 1:-1] + vn[0:-2, 1:-1])) + 
                    fy[1:-1, 1:-1] * dt)

    #conditions aux limites en xmin
    u[1:-1, 0] = (un[1:-1, 0] - un[1:-1, 0] * dt / dx *
                 (un[1:-1, 0] - un[1:-1, -1]) -
                 vn[1:-1, 0] * dt / dy * 
                 (un[1:-1, 0] - un[0:-2, 0]) - 
                 dt / (2 * rho * dx) * 
                 (p[1:-1, 1] - p[1:-1, -1]) + 
                 nu * (dt / dx**2 * 
                 (un[1:-1, 1] - 2 * un[1:-1, 0] + un[1:-1, -1]) +
                 dt / dy**2 *
                 (un[2:, 0] - 2 * un[1:-1, 0] + un[0:-2, 0])) + 
                 fx[1:-1, 0] * dt)
                
    v[1:-1, 0] = (vn[1:-1, 0] - un[1:-1, 0] * dt / dx *
                 (vn[1:-1, 0] - vn[1:-1, -1]) -
                 vn[1:-1, 0] * dt / dy *
                 (vn[1:-1, 0] - vn[0:-2, 0]) -
                 dt / (2 * rho * dy) * 
                 (p[2:, 0] - p[0:-2, 0]) +
                 nu * (dt / dx**2 * 
                 (vn[1:-1, 1] - 2 * vn[1:-1, 0] + vn[1:-1, -1]) +
                 dt / dy**2 * 
                 (vn[2:, 0] - 2 * vn[1:-1, 0] + vn[0:-2, 0])) + 
                 fy[1:-1, 0] * dt)

    #conditions aux limites en xmax
    u[1:-1, -1] = (un[1:-1, -1] - un[1:-1, -1] * dt / dx * 
                  (un[1:-1, -1] - un[1:-1, -2]) -
                  vn[1:-1, -1] * dt / dy * 
                  (un[1:-1, -1] - un[0:-2, -1]) -
                  dt / (2 * rho * dx) *
                  (p[1:-1, 0] - p[1:-1, -2]) + 
                  nu * (dt / dx**2 * 
                  (un[1:-1, 0] - 2 * un[1:-1,-1] + un[1:-1, -2]) +
                  dt / dy**2 * 
                  (un[2:, -1] - 2 * un[1:-1, -1] + un[0:-2, -1])) + 
                  fx[1:-1, -1] * dt)
                   
    v[1:-1, -1] = (vn[1:-1, -1] - un[1:-1, -1] * dt / dx *
                  (vn[1:-1, -1] - vn[1:-1, -2]) - 
                  vn[1:-1, -1] * dt / dy *
                  (vn[1:-1, -1] - vn[0:-2, -1]) -
                  dt / (2 * rho * dy) * 
                  (p[2:, -1] - p[0:-2, -1]) +
                  nu * (dt / dx**2 *
                  (vn[1:-1, 0] - 2 * vn[1:-1, -1] + vn[1:-1, -2]) +
                  dt / dy**2 *
                  (vn[2:, -1] - 2 * vn[1:-1, -1] + vn[0:-2, -1])) + 
                  fy[1:-1, -1] * dt)
    
    #on élimine les mauvaises valeurs à l'intérieur des surfaces
    for i in range(nx):
        for j in range(ny):
            if s[i,j]!=0:
                u[i,j]=0
                v[i,j]=0
    
    #vitesses initiales uniquement selon x
    u[:,0]=u0
    v[:,0]=0
    
    #test et actualisation de la boucle
    diff_u=(np.sum(u)-np.sum(un))/np.sum(u)
    stepcount=stepcount+1
            
#graphe
X,Y=np.meshgrid(x,y)
color=np.hypot(u,v) #méthode pour colorier les vecteurs
plt.quiver(X,Y,u,v,color)
plt.xlabel('X (m)')
plt.ylabel('Y (m)')
plt.title('['+stype.upper()+']'+'   Traînée: '+
          str(cx) + ' Dérive: '+str(cy)+' Compteur: '
          +str(stepcount),fontsize=10)
plt.suptitle('Représentation graphique des courants',fontweight='bold')
plt.colorbar(label='V (m/s)')
plt.scatter(X,Y,s,c='black',marker='s',alpha=0.5) #surface
plt.legend()
plt.savefig(f_name+'_'+stype+'_'+str(nx)+'.jpg',dpi=300)
plt.show()
print('process ended, graphic saved as '+f_name+'_'+stype+'_'+str(nx)+'.jpg')
Nous obtenons ainsi les figures suivantes:
Figure 6 - Cas du cube avec une discrétisation du plan par une matrice de taille (30×30)
Figure 7 - Cas du cube (100×100)
Figure 8 - Cas de la sphère (30×30)
Figure 9 - Cas de la sphère (100×100)
Conclusion: Ces résultats numériques sont en accord avec nos résultats théoriques et expérimentaux: les flux d'air contournent l'obstacle et les vitesses augmentent naturellement en évitant ces objets. Nous pouvons noter que l'effet “d'évitement” est plus marqué dans le cas du cube, ceci étant dû à une surface de référence plus brute. — === Seconde approche: programme(s) Java par des universitaires === Afin de comparer nos résultats, on se propose d'utiliser des logiciels existants afin de modéliser les lignes de courants de différents profils. Pour notre étude, nous avons choisis d'utiliser l'outil JavaFoil, soufflerie numérique 2D gratuite écrite et diffusée par le chercheur et modéliste allemand Martin HEPPERLE. Celle-ci permet de comparer rapidement différents profils en fonction des nombreux paramètres mis à notre disposition (nombres de Reynolds et de Mach, influence du sol, angles d'attaques, densité,…) et d'en tirer des informations notamment les coefficients de traînées, de portances et les fameuses lignes de courants (et bien d'autres, assez superficielles dans le cadre de notre étude). Voici un aperçu de l'interface du logiciel :
Figure 10 - Interface de Javafoil, vu de l'onglet Géométrie
Dans le cadre de notre étude, nous nous intéresserons principalement aux onglets “Géométrie”, “Modifications”, “Écoulements” et “Options”. Ainsi, passons à l'étude de l'aérodynamisme d'une sphère dans deux situations différentes: entourée d'un même fluide (ici de l'air) et à proximité du sol dans ce même fluide. On pose un nombre de Reynolds égal à $Re = 10^5$, un nombre de Mach égal à $Ma=0$ (l'objet d'étude est immobile) et enfin on considère un angle d'attaque égal à $\theta =0$ également (Remarque: l'angle d'attaque n'est pas vraiment un paramètre pertinent dans le cas d'une sphère, mais l'on doit bien choisir une valeur). On obtient les résultats suivants:
Figure 11 - Sphère entouré de fluide
Figure 12 - Considération des effets du sol pour cette même sphère
Légende: Les lignes blanches correspondent aux lignes “iso-Cp” c'est-à-dire les lignes où le coefficient de pression $C_p$ est constant.
  • *Observations Figure 11 : On obtient des résultats similaires à ceux obtenus avec notre première approche: notre démarche a donc été efficiente. De plus, JavaFoil nous indique que le coefficient de traînée pour ce profil est de $C_x = 0,471$ ce qui correspond aux valeurs données précédemment. Figure 12 : On constate des lignes de courants similaires avec néanmoins une turbulence beaucoup plus grande observée à proximité du sol et également une vitesse bien plus importante du fluide à proximité du sol. Afin de compléter notre étude des lignes de courants dans le cadre de la sphère, on se propose de représenter la distribution de coefficient de pression $C_p$ dans les mêmes conditions. Mais tout d'abord, présentons rapidement ce coefficient de pression. Il est définit par la formule ci-dessous: $$C_p = \frac{P- P_{air} }{\frac{1}{2} \rho_{air}V^2 }$$ avec $P$ la pression au point de mesure, $P_{air}$ la pression de référence (il s'agit souvent de la pression statique mesurée en amont) et $V$ la vitesse libre du fluide. Remarque relative à $C_p$: Le coefficient de pression ne peut jamais avoir une valeur supérieur à $1$ (limite fixé par le théorème de Bernoulli) par contre ce même coefficient ne possède, en théorie, pas de borne inférieure. Par ailleurs, la valeur de $C_p$ est essentiellement lié à la géométrie de la structure et reste relativement indépendante de la vitesse de l'écoulement. Représentons maintenant la distribution du coefficient de pression $C_p$ :
    Figure 13 - Distribution de $C_p$ pour une sphère entourée d'air
    On remarque ici une distribution inverse à celle de la vélocité: les zones de faibles vélocités correspondent aux zones ou la valeur de $C_p$ est maximale (environ $1$ ici) et à l'inverse, les zones de fortes vélocités correspondent aux zones où la valeur de $C_p$ est minimal (de l'ordre de $-2$) Enfin afin de diversifier notre étude, on se propose d'étudier un profil aérodynamique typique utilisé en aviation: un profil d'aile de type NACA (National Advisory Committee for Aeronautics). Ces profils possèdent des coefficients de traînées $C_x$ très bas et d'autres propriétés aérodynamiques utiles à l'aviation. À titre de comparaison, le $C_x$ d'une aile d'avion est compris entre $0,005$ et $0,010$ contre des valeurs comprises entre $0,5$ et $0,2$ dans le secteur de l'automobile ou encore pour une valeur de $1,2$ pour un corps humain debout. Ainsi pour notre étude, nous avons choisi le profil NACA 63-131. Son code nous donne les informations caractéristiques du modèle : il appartient à la série 6 de profil NACA, possède une pression minimum à la corde de 30 %, un coefficient de portance de $C_z=0,1$ et une épaisseur maximale de 31 %. On représente donc la distribution de vélocité pour ce profil, entouré d'air, via encore une fois notre outil favori JavaFoil.
    Figure 14 - Distribution de vélocité pour une aile d'avion entourée d'air
    On remarque ici un champ de vitesse beaucoup plus homogène qu'avec les objets étudiés précédemment: la formation de turbulence est limitée et des faibles augmentations de vitesses sont observés au-dessus et en-dessous du profil d'aile. On en déduit qu'une aile de ce modèle associe une faible traînée et une forte portance, deux caractéristiques optimales pour l'aviation, domaine pour lequel elle a été confectionnée. Enfin pour conclure cette étude via Javafoil, on se propose d'étudier l'écoulement pour une Formule 1. Le profil de celle-ci étant trop complexe pour le logiciel et impossible à retranscrire dans tous ces détails, on utilisera un simple triangle rectangle faisant office de Formule 1 “très simplifiée”. Évidemment pour cette étude, on prendra en compte les effets du sol correspondant ici à la bordure inférieur de l'interface (soit $y = 0$)
    Figure 15 - Distribution de vélocité pour une Formule 1 (approximée à un triangle) posée sur le sol
    On obtient ici un écoulement pour le moins intéressant: la vélocité est plus importante au devant de la Formule 1 (au niveau de l'aileron avant donc) et également au derrière de celle-ci, plus particulièrement à la place qu'occuperait l'aileron arrière. De ce fait, cette représentation souligne bien l'importance primordiale de l'aileron avant/arrière dans une Formule 1, deux variables donc où les constructeurs peuvent optimiser l'aérodynamisme de leurs véhicules. Enfin, on constate également une vélocité très importante sous le châssis de la voiture: ce phénomène souligne cette fois l'importance de la portance dans une Formule 1. Sans une portance élevé, l'adhérence, et donc la vitesse, se verrait hautement diminué pouvant allant jusqu'au décollage de la Formule 1 ! — === Bilan des deux méthodes === Au final, nous pouvons noter que le résultat obtenu avec notre code est plutôt similaire à celui obtenu par des logiciels plus poussés ou encore des simulations détaillées. Nous notons que plusieurs informations permettent de caractériser les lignes de courant comme l'intensité des vitesses en ces points et aussi différents coefficients comme celui de pression $C_p$. Cependant nos résultats obtenus à travers l'expérience sont loin de la théorie et ceci est du au manque de matériel précis mis à notre disposition et également les approximations que nous furent obligés de faire du coup. Ceci donne des coefficients de trainée qui ne correspondent pas aux théoriques même avec les incertitudes. On en déduit que notre méthode expérimentale n'est pas assez précise afin de mesurer des coefficients dépendant de tant de variables et paramètres. ===== Bibliographie =====     Articles scientifiques: - D'HONDT M., Étude théorique, expérimentale et numérique de l'écoulement de refroidissement et de ses effets sur l'aérodynamique automobile, Thèse de doctorat sous la direction de DEVINANT P., Orléans, École doctorale de sciences et technologies, 2010 [NNT 2010ORLE2026] - LECLERC C., Réduction de la traînée d'un véhicule simplifié à l'aide du contrôle actif par jet synthétique, Thèse de doctorat sous la direction de KOURTA A., Toulouse, Institut de Mécanique des Fluides de Toulouse, 2008 [NNT 2008INPT012H] Images externes (par ordre d'apparition): - CFD for Formula 1, http://www.4erevolution.com/wp-content/uploads/2015/09/Williams-F1-2.jpg [consulté le 22/02/17] - F1 Championnat, http://www.larevueautomobile.com/images/image-actu/F1-Resultat-1.jpg [consulté le 23/02/17] - Quelques profils et leurs traînées aérodynamiques, https://upload.wikimedia.org/wikipedia/commons/thumb/2/26/Drag-fr.svg/220px-Drag-fr.svg.png [consulté le 24/02/17] Liens Internet: - Notions d'aérodynamisme, site de LEDUC C., http://www.laberezina.com/technique/aerodynamique.htm [en ligne, consulté le 11/02/17] - Aérodynamisme, page de WIKIPEDIA, https://fr.wikipedia.org/wiki/A%C3%A9rodynamique [en ligne, consulté le 24/02/17] - Aérodynamisme automobile, page de WIKIPEDIA, https://fr.wikipedia.org/wiki/A%C3%A9rodynamique_automobile [en ligne, consulté le 11/02/17] - Coefficient de traînée, page de WIKIPEDIA, https://fr.wikipedia.org/wiki/Coefficient_de_tra%C3%AEn%C3%A9e [en ligne, consulté le 24/02/17] - Aérodynamisme d'une Formule 1, page de MADIER D., http://www.f1-forecast.com/pdf/F1-Files/F1-Forecast%20-%20Aerodynamique%20Formule%201.pdf [en ligne, consulté le 28/02/17] - CFD Python - 12 steps to Navier-Stokes, page de BARBA L. (mathématicienne, professeur à l'Université George Washington), http://lorenabarba.com/blog/cfd-python-12-steps-to-navier-stokes/ [en ligne, consulté le 01/03/17] - Ahmed Body, page de CFD Online, https://www.cfd-online.com/Wiki/Ahmed_body [en ligne, consulté le 05/03/17] - Notions d'aérodynamique, cours de l'École Polytechnique, http://www.editions.polytechnique.fr/files/pdf/EXT_1332_5.pdf [en ligne, consulté le 05/03/17] Ouvrages**:
  1. AMIROUDINE S., BATTAGLIA J.-L., Mécanique des fluides, Paris, éd. Dunod, 2011 [ISBN 978-2-10-054933-7]
  2. FAURE T., Dynamique des fluides appliquée : applications à l'aérodynamique, Paris, éd. Dunod, 2008, coll. “Sciences Sup” [ISBN 978-2-10-051099-3]
  3. FLETCHER C.A.J., Computational Techniques for Fluid Dynamics, Berlin, éd. Springer-Verlag, 1991, 2e édition [ISBN vol. I 3-540-53058-4, ISBN vol. II 3-540-53601-9]
  4. GUYON P., HULIN J.-P., PETIT L. et al., Physical Hydrodynamics, Oxford, éd. Oxford University Press, 2001 [ISBN 0-19-851745-9]
  5. HUCHOT W.-H., Aerodynamics of road vehicles, éd. Society of Automotive Engineers, 1998, 4e édition [ISBN 0-7680-0029-7]
  6. OUZIAUX R., PERRIER J., “Dynamique des fluides compressibles” (ch. 5) dans Mécanique des fluides appliquée, Paris, éd. Dunod, 1998, 3e édition, coll. “Sciences Sup” [ISBN 2-10-048400-1]
  7. POPE S.B., Turbulent Flows, Cambridge, éd. Cambridge University Press, 2010, 7e édition [ISBN 978-0-521-59886-6]
  8. RIEUTORD M., Une introduction à la dynamique des fluides, Bruxelles, éd. De Boeck, 2014 [ISBN 978-2-8041-8554-1]