-- Simulation numérique en mécanique des fluides --

TP n°6

 Méthode des singularités

(méthode des panneaux)
simon.marie@lecnam.net
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
# la librairie time nous sera utile pour calculer le temps d'execution des scripts
import time 
#Option d'affichage et taille de police sur les figures:
fs=20
plt.rc('xtick',labelsize=fs)
plt.rc('ytick',labelsize=fs)
%matplotlib inline

On se propose dans ce TP, de se familiariser avec la méthode des singularités (ou méthode des paneaux). On utilisera pour ce TP le profil NACA0012 mais vous pourrez changez le profil une fois les bon résultats obtenus. Vous pouvez utiliser pour cela les données du site Airfoiltools.

La méthode des panneaux consiste à discrétiser une surface en un ensemble discret de panneaux puis de décrire chaque panneaux à l'aide d'un potentiel complexe. Le potentiel général prenant en compte un écoulement uniforme et l'influence de chaque paneau est alors écrit sous la forme: $$ \Phi(x,y)=U_0(x\cos(\alpha))+y\sin(\alpha))+\sum_{j=1}^{N}\int_{L_j}\left[\dfrac{q_j(s)}{2\pi}\ln(r(s))-\dfrac{\gamma}{2\pi}\theta_j(s)\right]ds $$ Avec $N$ le nombre de panneaux, $r(s)=\sqrt{[x-x_j(s)]^2+[y-y_j(s)]^2}$

et $\theta_j(s)=arctan\left(\dfrac{y-y_j(s)}{x-x_j(s)}\right)$

Cette méthode modélise donc chaque panneau à l'aide d'un ensemble de sources et de tourbillons. Dans l'approche utilisée ici, les sources ont des intensité $q_j$ différente sur chaque panneau alors que l'intensité du tourbillon $\gamma$ est la meme sur tous les panneau.

Les inconnues du problème sont donc les $N$ coefficients $q_j$ (les sources des chaques panneaux) et le coefficient $\gamma$ (un vortex identique sur chaque panneau). Pour calculer ces $N+1$ inconnues, on doit trouver un système linéaire à $N+1$ équations grâce aux conditions suivante:

  • Condition 1: Vitesse tangentiel nulle sur chaque panneaux: $\mathbf{U_i(x_c,y_c)}.\mathbf{n_i}=0$ (N équations)
  • Condition 2: Pression identique au bord de fuite (Kutta): $\mathbf{U_1(x_c,y_c)}.\mathbf{t_1}=-\mathbf{U_N(x_c,y_c)}.\mathbf{t_N}$ (1 équation)

Avec $\mathbf{n_i}=\left(\begin{array}[l] --\sin(\theta_i)\\ \cos(\theta_i)\end{array}\right)$ le vecteur normal à chaque panneaux et $\mathbf{t_i}=\left(\begin{array}[l] \\cos(\theta_i)\\ \sin(\theta_i)\end{array}\right)$ le vecteur tangentiel à chaque panneau.

Travail préliminaire

Afin de réaliser ce TP, vous devez avoir correctement terminé l'exercice des TD consacré aux écoulements potentiels et à la méthode des panneaux.

1 - Maillage de la surface

La première étape consiste à choisir une surface à étudier. Ici en 2D on choisir le profil d'aile NACA0012 qui est un profil symétrique.

On utilisera pour cela les données fournies dans le fichier NACA012.dat. On prendra dans un premier temps $U_0=1.$ et $\alpha=0$ (incidence nulle) puis une fois le code validé, on prendra $\alpha=10°$.

In [ ]:
U0=1.
alpha=0.*np.pi/180.

Charger les données du fichier NACA012.dat:

In [ ]:
# Téléchargement des donnes pour la première execution:
#! wget https://hpp.education/Lessons/CFD_FIP/Files/NACA012.dat
In [ ]:
X=np.loadtxt('NACA012.dat')[::-1]
x,y=X[:,0],X[:,1]

Fixer le nombre de panneau en se basant sur le nombre de point décrivant le profil dans le fichier dat:

In [ ]:
N=... # Choisir le nombre de panneaux
print(N)

On pourra représenter le profil en soulignant le premier point et le dernier point.

In [ ]:
fig=plt.figure(figsize=(16,4))
plt.plot(x,y,'or')

# Premier point
plt.plot(x[0],y[0],'sg')

# Dernier point
plt.plot(x[-1],y[-1],'.k')

plt.axis('equal');plt.grid(True)

#Zoom sur le bord de fuite:
#plt.xlim(0.95,1.05);plt.ylim(-0.01,0.01)

Commentaires:

-

2 - Calcul des panneaux

La deuxième étape consiste à calculer les coordonnées des panneaux:

  1. Longueur $L_i=\sqrt{ (x_{i+1}-x_i)^2+(y_{i+1}-y_i)^2 }$

  2. Position des centres ($xc_i,yc_i$)

  3. Orientation ($\theta_i$). On utilisera la fonction arctan2 qui prend bien deux argument en entrée. (Cf doc)

In [ ]:
# Calcul de la longueur de chaque paneaux et des points centraux:
L,xc,yc,theta2=np.zeros(N),np.zeros(N),np.zeros(N),np.zeros(N)

# Distance entre les bords de deux panneau:
dpx=...
dpy=...

# Longueur du paneau
L[:]=np.sqrt(dpx**2+dpy**2)

# Centre du panneau
xc[:]=(x[1:]+x[0:-1])/2
yc[:]=(y[1:]+y[0:-1])/2

# Orientation
theta2[:]=np.arctan2(dpy,dpx)

# On représente theta en fonction de la longueur des panneau:
fig=plt.figure(figsize=(16,4))
plt.plot(L,theta2,'o')

Commentaires:

-

3 - Calcul des coefficients de la matrice $\mathbf{A}$

La troisième étape consiste à calculer les coefficients de la matrice $A$ et du vecteur $B$. Ces coefficients viennent des conditions aux limites imposées (vitesse tangentielle et contition de Kutta) aux points centraux des panneaux.

A partir des calculs réalisés en TD, on peut écrire les 2 conditions sous la forme d'un système linéaire $\mathbf{A.Q}=\mathbf{B}$. Aevc:

$$ \mathbf{A}= \left[ \begin{array}{cccc} a_{11}&\cdots&a_{1N}&\sum_j b_{1j}\\ \vdots&\ddots&\vdots&\vdots\\ a_{N1}&\cdots&a_{NN}&\sum_j b_{Nj}\\ -b_{11}-b_{N1}&\cdots&-b_{1N}-b_{NN}&\sum_j a_{1j}+a_{Nj} \end{array} \right] $$

et

$$ \mathbf{B}=[\cdots,-U_0\sin(\alpha-\theta_i),\cdots,-U_0(\cos(\alpha-\theta_0)+\cos(\alpha-\theta_N))]^T $$

avec:

$$ a_{ij}=\int_{0}^{L_j}\dfrac{[yc_i-y_j(s)]\cos(\theta_i) - [xc_i-x_j(s)]\sin(\theta_i)}{[xc_i-x_j(s)]^2+[yc_i-y_j(s)]^2}ds $$$$ b_{ij}=\int_{0}^{L_j}\dfrac{-[yc_i-y_j(s)]\sin(\theta_i) - [xc_i-x_j(s)]\cos(\theta_i)}{[xc_i-x_j(s)]^2+[yc_i-y_j(s)]^2}ds $$

Pour calculer les coefficients de la matrice, on peut définir une fonction générique:

$$ I_{ij}(X,Y)=\dfrac{1}{2\pi}\int_{0}^{L_j}\dfrac{[xc_i-x_j(s)]X - [yc_i-y_j(s)]Y}{[xc_i-x_j(s)]^2+[yc_i-y_j(s)]^2}ds $$

que l'on pourra calculer à l'aide de la fonction quad du module integrate.

In [ ]:
# On va définir une fonction qui calcul l'intégrale pour nous.
# Cette fonction va nous servir à la fois pour les sources et pour les vortex.

def integral(xi,yi,xj,yj,thetaj,Lj,dx,dy):
    def integrand(s):
        xs=...
        ys=...
        XY=(xi-xs)**2+(yi-ys)**2
        return ((yi-ys)*dy+(xi-xs)*dx)/XY
    return integrate.quad(integrand,0.,Lj)[0]
In [ ]:
Aij=np.zeros((N+1,N+1)) # Matrice totale du système
a=0.5*np.eye(N) # Matrice des sources => qj
b=np.zeros((N,N)) # Matrice des vortex  => gamma

Bi=np.zeros(N+1) # Vecteur de la partie droite du système.

for i in np.arange(N):
    for j in np.arange(N):
        if i!=j:
            a[i,j]=integral(???)/(2.*np.pi)
            b[i,j]=integral(???)/(2.*np.pi)
    Bi[i]=???

Aij[:N,:N]=???
Aij[:N,-1]=???
Aij[-1,:N]=???
Aij[-1,-1]=???
Bi[-1]=???

Commentaires:

-

4 - Résolution du système linéaire

La quatrième étape consiste à résoudre le système linéaire:

$$ \mathbf{A}Q=\mathbf{b} $$

On utilisera pour cela la fonction solve de numpy.linalg.

In [ ]:
Q=np.linalg.solve(Aij,Bi)

On peut vérifier la précision de la solution en calculant la somme des sources qui doit être nulle: $$ \sum_i^N q_jL_j=0 $$

In [ ]:
plt.plot(Q[:N])
print(np.sum(Q[:N]*L))

Commentaires:

5 - Calcul de la vitesse tangentielle

A partir des solutions, calculer la vitesse tangentiel obtenue sur chacun des points centraux des panneaux.

La vitesse tangentiel sur le paneau $i$ s'écrit:

$$ U_{t,i}=U_x(x_{c,i},y_{c,i})\cos(\theta_i)+U_y(x_{c,i},y_{c,i})\sin(\theta_i) $$

En remplacant on obtient:

$$ U_{t,i}=U_0\cos(\theta_i-\alpha) + \sum_j I_{ij}(\cos(\theta_i),\sin(\theta_i))q_j + \sum_j I_{ij}(-\sin(\theta_i),\cos(\theta_i))\gamma $$

Avec $I_{ij}(X,Y)$, la fonction intégrale définie précédemment:

In [ ]:
Ut=np.zeros(N)
for i in np.arange(N):
    Ut[i]=U0*np.cos(theta2[i]-alpha)
    for j in np.arange(N):
        if i!=j:
            Ut[i]+=...
            Ut[i]+=...
In [ ]:
plt.plot(xc,Ut)

Commentaires:

-

6 - Calcul des coefficients de pression $C_p$ et de portance $C_L$

A partir des résultats, on peut tracer le coefficient de pression:</h1>

Enfin à partir des résultats, on peut tracer le coefficient de pression: $$ C_p=1-\dfrac{U_t^2}{U_0^2} $$

Tracer cette évolution en fonction de la position sur le profil et compare aux valeurs de références fournies dans les fichiers suivants:

In [ ]:
# Téléchargement des donnes pour la première execution:
#! wget https://hpp.education/Lessons/CFD_FIP/Files/Cp_NACA12_0deg.dat
#! wget https://hpp.education/Lessons/CFD_FIP/Files/Cp_NACA12_10deg.dat
In [ ]:
# Chargement du Cp de référence
Cpt=np.loadtxt('../Cp_NACA12_0deg.dat')
#Cpt=np.loadtxt('../Cp_NACA12_10deg.dat')

#Calcul du Cp à partir de la vitesse:
Cp=...

fig=plt.figure(figsize=(16,4))
plt.plot(xc,Cp,'-o')
plt.plot(Cpt[:,0],Cpt[:,1],'rs')
plt.grid(True)

Commentaires:

-

La condition de Kutta-Joukowski permet d'écrire le coefficient de portance: $$ C_L=\dfrac{2\Gamma}{U_0 c} $$ Ou $\Gamma$ est la circulation totale donnée par $\displaystyle \Gamma=\sum_{i=1}^{N}\gamma L_i$

Calculer le coefficient de portance ainsi obtenu:

In [ ]:
Cl=...
print(Cl)

Commentaires:

-

7 - Maillage du plan et ligne de courant

On peut maintenant représenter l'écoulement autour du profil en évaluant les vitesses $U_x$ et $U_y$ sur les points d'un maillage régulier. On prendra comme domaine $(x,y)\in [-2,3],[-1,1]$ avec un pas de grille de $\Delta x=0.1$.

In [ ]:
X=np.arange(-2,3,0.1)
Y=np.arange(-1,1,0.1)
A=np.zeros((len(X),len(Y),N))
B=np.zeros((len(X),len(Y),N))
Atx=np.zeros((len(X),len(Y),N+1))
Aty=np.zeros((len(X),len(Y),N+1))

for i, xi in enumerate(X):
    for j,yj in enumerate(Y):
        for k in np.arange(N):
            if i!=j:
                A[i,j,k]=...
                B[i,j,k]=...


Atx[:,:,:-1]=A
Atx[:,:,-1]=np.sum(B,axis=2)
−
Aty[:,:,:-1]=B
Aty[:,:,-1]=np.sum(A,axis=2)

bx=U0*np.cos(alpha)*np.ones((len(X),len(Y)))
by=U0*np.sin(alpha)*np.ones((len(X),len(Y)))
In [ ]:
Ux=np.dot(Atx,Q)+bx
Uy=np.dot(Aty,Q)+by
In [ ]:
fig=plt.figure(figsize=(15,6))
plt.pcolormesh(X,Y,Ux.T,shading='gouraud',cmap='jet')
plt.colorbar()
Xx,Yy=np.meshgrid(X,Y)
plt.streamplot(Xx,Yy,Ux.T,Uy.T,2)
plt.plot(x,y,'k')

On pourra reprendre les étapes du TP en changeant le nombre de panneaux, l'incidence ou même prendre un profil différent.

Conclusion

Présenter une synthèse de vos résultats et de votre démarche.

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

Retour en haut de la page