INTRODUCTION

Depuis Otto Lilienthal et Clément Adler dans les années 1890, les engins volants on énormément évolué.

C’est particulièrement vrai du point de vue des performances, car entre les 50 mètres de Clément Adler, et les milliers de kilomètres que font tous les jours nos longs courriers, Il a fallu modifier et optimiser énormément de caractéristiques. L’une d’entre elles, parmi les plus importantes, c’est l’aérodynamique.

Chacun a vu dans son enfance, que le fait de passer la main par les vitres d’une voiture engendrait une résistance. L’air en mouvement apporte une pression sur la main. Le vol est basé sur ce principe : s’appuyer sur l’air. Tout l’art, maintenant, étant de maximiser ce facteur d’appui : la portance, tout en réduisant la résistance de l’air au vol qui engendre des pertes plus que substantielles d’énergies, et donc de carburant.

Il est nécessaire pour cela d’étudier la forme des ailes : que se soit la forme projetée dans le plan : delta, trapèze, …, et surtout, ce qui est particulièrement important : le profil (section) de l‘aile.

Divers moyens, expérimentaux, ou calculatoires existent pour ces études. L’essor de l’informatique de ces 30 dernières années permet à moindre coût de faire des études aérodynamiques numériques fiables, pour autant que l’on maîtrise les hypothèses utilisées dans le code de calcul.

Une des méthodes utilisée est « la méthode des panneaux ». Elle a servi au développement des avions commerciaux actuels, dont le prestigieux concorde, unique avion commercial supersonique.

Or, cette méthode qui est encore à la base de codes commerciaux utilisés en construction aéronautique existe sous plusieurs formes ayant chacune avantages et inconvénients. L’étude qui va suivre permettra de faire le point sur toutes ces variantes utilisées dans l’étude de profil d’ailes, et développera une de ces méthodes : la méthode de panneaux avec tourbillons par formulation de ligne de courant.

Les résultats mettront notamment en évidence l’importance du nombre de panneaux dans la précision du calcul, les répartitions des coefficients de pression sur le profil, ainsi que l’évolution de la portance en fonction de l’incidence.

La pression et la vitesse sur un profil d’aile, sont obtenues par la résolution des équations de Navier Stokes. Or, étant donné que la résolution de ces équations est très compliquée pour un fluide réel, on est obligé de poser des hypothèses pour simplifier le problème. On va donc considérer que le fluide qui s’écoule autour du profil, est un fluide parfait, donc non visqueux et irrotationnel. Ces hypothèses sont bien sûr justifiables ! En effet, comme la vitesse d’écoulement est rapide (bien que subsonique), l’effet des forces de viscosité est négligeable par rapport aux forces d’inertie. De plus, l’écoulement est irrotationnel si la vorticité est nulle en tout point. (Ñ ´ = 0) Et enfin, puisque le fluide est incompressible, la densité est constante tout au long de sa trajectoire. Cela nous permet d’obtenir des simplifications des équations de continuité et du mouvement. De telle sorte que :

D’autre part, étant donné que l’écoulement est irrotationnel, on sait que = ÑF et on en déduit que l’équation de Laplace est de la forme :

Les solutions de cette équation fournissent le potentiel de vitesse d’un écoulement irrotationnel et incompressible. Puis l’équation du mouvement nous permet de trouver la pression. L’équation de Bernoulli s’obtient de l’intégrale de l’équation du mouvement, en tenant compte que la densité est constante. Ainsi :

En utilisant l’équation de Bernoulli, on obtient le coefficient de pression suivant :

Le niveau de pression de référence n’influe pas directement sur la traînée car la force de pression sur le corps étudié est nulle seulement si la pression est constante sur la totalité de la surface étudiée.

L’équation de Laplace nous permet donc de résoudre l’écoulement irrotationnel et incompressible autour de l’aile de l’avion. Ainsi, la détermination d’une fonction harmonique F doit satisfaire les conditions aux limites sur les frontières du domaine pour répondre au problème d’écoulement.

Nous admettons que les conditions aux limites correspondent à l’infini, où la perturbation occasionnée par le profil, tend vers zéro, et sur la surface du corps S_c, où l’on néglige la condition d’adhérence puisque le fluide est considéré comme parfait, donc non visqueux. De plus on suppose que la surface S_c est imperméable, donc que l’écoulement n’y pénètre jamais. Toutes ces hypothèses étant posées, on va maintenant tenter de résoudre le problème mathématique suivant :

Les équations (7) et (8) sont les conditions frontières qui seront utilisées lors de la résolution de l’équation de Laplace. Notons que Ñ²Y = 0 satisfait aussi l’équation de Laplace.

En théorie, on peut avoir une infinité de solutions pour l’intensité de la circulation, de sorte que le potentiel total de l’écoulement satisfait aux conditions de continuité, d’irrotationnalité et de l’écoulement tangent pour un profil dans un écoulement stationnaire à un angle d’attaque donné. On comprend donc que ces conditions sont insuffisantes pour trouver une solution unique du problème. Il est donc nécessaire de poser une condition supplémentaire pour obtenir une valeur unique de la circulation. Cette condition est issue de l’étude de l’écoulement au bord de fuite réalisée par Wilhelm Kutta et est formulée de la manière suivante :

· Pour un profil se terminant par un point de rebroussement, la vitesse à l’extrados et à l’intrados du bord de fuite, a la même valeur finie. Dans ce cas, les deux vitesses ont même direction et même valeur finie. Cependant, la pression au bord de fuite a une valeur unique p_BF et l’équation de Bernoulli appliquée à l’extrados et à l’intrados du profil à proximité du bord de fuite, donne :

· Pour un profil se terminant par un dièdre, le bord de fuite est un point d’arrêt. Soit les vitesses de l’écoulement au bord de fuite, v_e tangente à l’extrados et v_i tangente à l’intrados ; on a deux vitesses de directions différentes en un même point, ce qui n’est possible que si les deux vitesses sont nulles et que le point en question est un point d’arrêt. D’où v_e = v_i = 0

En résumé, les conditions aux frontières permettent de résoudre l’équation de Laplace, le champ de vitesse est déterminé à l’aide de l’équation qui montre que la vitesse dérive d’un potentiel. La distribution de pression est calculée avec l’équation (5) du coefficient de pression. Le problème d’écoulement irrotationnel et incompressible est résolu et, puisque l’équation de Laplace est une équation linéaire aux dérivées partielles du second ordre, on peut appliquer le principe de superposition des solutions. L’écoulement est représenté par la somme des écoulements incompressible irrotationnels élémentaires.

Les méthodes des panneaux sont des méthodes d’étude numérique des profils d’ailes. Il en existe plusieurs:

Le principe de discrétisation consiste à “facétiser” un profil, c’est à dire que l’on va le représenter par une série de petits segments de longueur variable appelés panneaux. Cette discrétisation ne doit pas être aléatoire et doit permettre de construire un polygone le plus proche possible du profil. De plus, il est nécessaire d’avoir une densité de panneaux plus importante aux extrémités du profil (bords d’attaque et de fuite) car ils constituent les endroits où la courbure est la plus prononcée, donc plus difficile à représenter par des segments. On utilise en général la formule suivante pour déterminer les coordonnées des panneaux :

x_j représente la distance du bord d’attaque au point j sur la corde et N, le nombre de panneaux. On prend comme point de contrôle, le point situé au centre du panneau. Les coordonnées de ce point sur le panneau i sont :

En général, le choix d’une soixantaine de panneaux fournit une bonne précision de la discrétisation. Bien sûr, ceci peut varier selon la cambrure des profils.

Après avoir réalisé la discrétisation, on doit procéder à la distribution des singularités sur les panneaux pour simuler l’écoulement du fluide. Cette répartition et les hypothèses de base associées peuvent varier d’une méthode à l’autre. Le principe reste qu’il faut poser des conditions qui limitent le nombre d’inconnues tout en fournissant des résultats valables et en simplifiant les calculs. Puis, le calcul de vitesse sur chaque panneau permet de trouver le coefficient de pression en chaque point duquel découle le coefficient de portance et la circulation.

Toutes les méthodes de panneaux reposent sur l’identité de Green. Celle ci permet de calculer la valeur du potentiel en tout point du champ concerné, pour autant que l’on connaisse les conditions aux limites. La surface en question, ne doit cependant contenir aucune singularité. Puisque dans notre cas, les singularités du profil sont distribuées à la surface du profil, examinons le champ décrit à la figure suivante.

Figure 4 Surface représentative à l’intérieur de laquelle l’identité de Green est applicable

La première intégrale de l’équation (13) représente l’influence au point P d’une distribution de sources sur la surface S et le terme est la vitesse normale à la paroi. Dans notre cas elle est nulle pour satisfaire la condition de tangence des vitesses. Donc le premier terme est nul. Dans la seconde partie de l’équation, le terme représente le taux de changement du potentiel F_sen direction . On peut montrer que ce terme correspond à l’influence d’un doublet.

Si on pose F₁ et F₂ les valeurs du potentiel aux points Q₁ et Q₂, on a alors :

Si on interprète F₁ et F₂ comme étant respectivement une source et un puits d’intensité unitaire, alors on obtient bien la définition du doublet. Selon Moran, la contribution de la surface S_¥ au point P(x_p, y_p) est donnée par :

En tenant compte de ce dernier terme, on peut réécrire l’identité de Green de la façon suivante :

Il existe une deuxième façon de considérer le problème. On pose que V = V_¥ + ÑF, où F est le potentiel de perturbation. La surface S_¥ ne contribue pas au potentiel, on a alors :

Le NACA (National Advisory Committee for Aeronautics), prédécesseur de la NASA, a étudié différentes familles de profiles répondant à diverses applications. Parmi celles ci, on peut distinguer la famille de profils à quatre chiffres, celle à cinq chiffres et les profils laminaires portant la désignation NLF. Dans ce rapport, nous allons principalement nous intéresser à ceux à quatre chiffres.

Dans cette famille, un profil est représenté par quatre chiffres. Le premier indique la cambrure maximale relative en pourcentage de la corde ; le deuxième représente la position de cette cambrure en pourcentage de la corde et les deux derniers spécifient l’épaisseur relative maximale en pourcentage de la corde. Par exemple, dans le cas du profil 4412 :

Nous allons maintenant voir les formules permettant de calculer les coordonnées des points d’un profil. Dans le cas des profils symétriques de cette famille, on se sert de l’expression suivante :

Pour le cas des profils cambrés, la forme géométrique de la cambrure moyenne est représentée par deux paraboles :

Où le point A (x_A, y_A), situé à la distance où la flèche est maximale, représente le point commun des deux paraboles. Remarquons que y_A représente la cambrure maximale et que x_A en donne la position.

Dans cette famille, cinq chiffres caractérisent les profils. Le premier représente le coefficient de portance caractéristique, les deux suivants indiquent la position de la cambrure maximale et les deux derniers spécifient l’épaisseur relative maximale.

· L’indice 2 est la marge au-dessus et au-dessous du coefficient de portance caractéristique pour laquelle il existe un gradient de pression favorable (0,2)

Cette méthode élaborée par Hess et Smith en 1966 consiste tout d’abord à diviser le profil en un certain nombre de panneaux afin d’en approcher le mieux possible le contour réel. On modélise ensuite l’écoulement à l’aide de sources et de tourbillons distribués sur chaque panneau. Ensuite, les conditions de tangences et de Kutta, nous permettent de déterminer l’intensité de ces singularités. Puis, une fois celles ci connues, on peut calculer la vitesse, et, par conséquent, le coefficient de pression en tout point de la surface, ainsi que la distribution de l’angle d’attaque en fonction de l’angle d’incidence du profil.

Pour un écoulement irrotationnel, le potentiel total autour du profil correspond à l’équation suivante :

Hess et Smith ont posé les conditions suivantes : l’intensité des sources est constante sur chacun des panneaux mais varie d’un panneau à l’autre, tandis que l’intensité des tourbillons est la même sur tous les panneaux. Donc sur un panneau i, q(s) = q_i tandis que g(s) = g. La discrétisation se fait de la façon suivante :

Avant d’effectuer les calculs, on doit préciser les conditions aux frontières. Tout d’abord, on considère que les points de contrôle sont situés au centre des panneaux. Puis on impose la condition de tangence entre les panneaux, ce qui donne une vitesse normale à la surface des panneaux nulle. Ensuite, la condition de Kutta établit que les vitesses tangentes du premier et du dernier panneau sont égales. L’ensemble de ces conditions permet de déterminer les inconnues de départ, c’est à dire l’intensité des sources et des tourbillons. Une fois le champ de vitesse connu, on peut en déduire le coefficient de pression sur chaque panneau. Puis le coefficient de portance est calculé à l’aide de la relation de Kutta Joukowski.

La méthode de Hess et Smith possède l’inconvénient de demander des calculs assez complexes et importants, notamment pour la condition de Kutta. Cependant c’est une méthode qui fournit de bons résultats en 2 dimensions comme en 3 dimensions.

Cette méthode vient de Mavriplis qui l’a présentée en 1971. Elle est basée sur le principe de superposition de solutions élémentaires consistant à distribuer des tourbillons sur la surface à étudier. Comme avec la méthode de Hess et Smith, on discrétise le profil à l’aide d’un ensemble de lignes droites que l’on nomme panneaux. On détermine la solution en distribuant des tourbillons sur la surface du corps plongé dans un écoulement uniforme avec un écoulement circulatoire. Pour ce type de représentation, Martensen a démontré que la condition de vitesse tangentielle nulle sur la surface interne du corps est équivalente à la condition de tangence et que la vitesse tangentielle sur la surface externe du corps est égale à l’intensité des tourbillons. A partir de ces observations, on peut exprimer la vitesse tangentielle en un point intérieur de la surface du corps, de la manière suivante :

L’équation (29) est nulle car le point P appartient à la surface du corps, donc la vitesse tangentielle est nulle.

On cherche donc la distribution g et l’écoulement circulatoire w/2 qui satisfont à l’équation (29) pour tous les points de la surface du corps. Une fois qu’on a déterminé la distribution g, on sait que la vitesse tangentielle en un point P de la surface externe du corps est égale à l’intensité du tourbillon à cet endroit.

Nous avons donc N inconnues g_j et une inconnue w. Les conditions de tangence nous fournissent N+1 inconnues et N équations. On va donc utiliser la condition de Kutta au bord de fuite comme dernière équation nécessaire à la résolution du problème. Comme la vitesse tangentielle en chaque point est donnée par l’intensité du tourbillon en cet endroit, la condition de Kutta se résume donc à :

Le système d’équations donne la valeur des vitesses aux points de contrôle, ce qui permet d’en déduire les coefficients de pression et de portance.

L’avantage de cette méthode est qu’elle est beaucoup plus facile à mettre en œuvre que la précédente.

Cette méthode se base également sur le principe de superposition de solutions élémentaires et utilise aussi des tourbillons. La vitesse tangentielle en un point est aussi donnée par l’intensité du tourbillon en ce point. La différence avec la méthode de formulation de vitesse réside dans l’équation intégrale qui n’exprime pas la vitesse tangentielle, mais plutôt la ligne de courant. La vitesse tangentielle sur la surface interne est nulle. La formulation de l’équation intégrale est :

Cette équation donne la fonction de courant en un point P quelconque. Puisque la fonction de courant est constante sur la surface du corps, pour tout point P de la surface du profil, on a donc :

En tenant compte de cette égalité et en remplaçant les fonctions de courant Y et Y_¥ par leurs valeurs respectives, la fonction intégrale sur le profil discrétisé en N panneaux, s’exprime comme suit :

On obtient N+1 inconnues : N valeurs de g_j et la valeur de Y_B. Or on a seulement N équations : la condition de tangence sur les N panneaux. La condition de Kutta au bord de fuite va nous fournir l’équation manquante pour solutionner le problème. Celle ci s’exprime de la même manière que dans le cas de la méthode des tourbillons avec formulation de vitesse, c’est à dire :

Après avoir calculé la vitesse tangentielle sur chaque point de contrôle (égale à l’intensité du tourbillon au même point), on peut en déduire les coefficients de portance et de pression.

Les principaux avantages de cette méthode sont qu’elle est assez simple et rapide. Mais elle ne peut s’utiliser que pour l’étude d’écoulements bidimensionnels.

Cette méthode utilise la théorie des fonctions de Green pour obtenir l’équation intégrale. Elle résulte du travail de Morino et Kuo et se base sur la solution de l’équation :

Elle consiste à supposer le potentiel constant sur chaque panneau, mais variant d’un panneau à l’autre. Le potentiel correspond à l’intensité des doublets. On pose que la circulation vaut :

On obtient donc un système de N inconnues, soit les potentiels F_j sur tous les panneaux. On peut remarquer que l’on a une inconnue de moins que dans les méthodes précédentes. On n’a donc besoin que de N équations, fournies par l’application de l’équation (32) au centre de chaque panneau. On obtient ainsi un système de N équations linéaires et N inconnues. La condition de Kutta est implicitement vérifiée. Pour obtenir les vitesses tangentielles, on dérive les potentiels. On obtient donc :

avec d, la distance entre le milieu des deux panneaux adjacents. Cette équation nous fournit la valeur de la vitesse tangentielle à chaque nœud du profil. Il ne nous reste plus qu’à calculer les coefficients de pression et de portance.

On peut obtenir la solution discrète du problème envisagé au moyen d’une méthode de panneaux développée par Moran, basée sur la solution de l’équation (30). La grande différence avec la méthode précédente, se situe sur le plan de la discrétisation de la distribution de doublets. En effet, Moran utilise une distribution linéaire de doublets. La distribution du potentiel sur un panneau devient donc :

j et x varient entre 0 et lj. La circulation est donc égale à la différence de potentiel entre le N+1^ème nœud et le 1^er nœud :

On applique maintenant l’équation de Green à un profil divisé en N panneaux. On obtient alors :

La résolution de l’intégrale permet de faire tendre le point P vers le i^ème nœud du profil. On obtient alors N+1 inconnues : les F_i, avec i = 1,…,N+1. On a donc un système d’équations algébriques. Mais, si on l’intègre, on fait apparaître des termes logarithmiques et si on applique juste l’équation (36) au point i, on obtient des singularités logarithmiques introduites dans l’équation. Il faut donc faire tendre le point P vers le i^ème nœud et observer ses effets sur l’équation de départ. On aboutit alors à :

Cependant, on peut calculer leur somme, qui est égale à l’angle entre les panneaux i-1 et i. Les termes réguliers sont tous les termes où j ¹ i et j ¹ i-1. Lorsque i = 1, on obtient un résultat particulier :

L’équation (36) fournit le même résultat lorsqu’elle est appliquée aux nœuds 1 et N+1. On a donc N équations pour N+1 inconnues. La condition de Kutta au bord de fuite va donc nous permettre d’obtenir un système de N+1 équations avec N+1 inconnues. Pour cela, on égalise les vitesses tangentielles au milieu des panneaux 1 et N, d’où :

La dérivée de l’équation de la distribution de potentiel sur le panneau permet de calculer la vitesse tangentielle au centre du panneau. On a donc :

Il suffit maintenant d’en déduire les coefficients de portance et de pression.

Cette méthode donne de bons résultats du fait que les points de contrôle soient situés aux nœuds et non pas en leur centre, mais les calculs qu’elle demande sont lourds.

Elle consiste en une division du profil en petits segments, les panneaux. La densité de panneaux doit être plus importante aux bords de fuite et d’attaque pour avoir une plus grande précision. Chaque panneau contient un tourbillon distribué uniformément mais dont l’intensité varie d’un panneau à l’autre.

La discrétisation en N panneaux entraîne N+1 inconnues, à savoir les N valeurs de gi et celle de YB. Après évaluation de l’intégrale, la fonction de courant devient :

Lorsque le point P est situé au centre du panneau j, il faut utiliser l’expression suivante :

On obtient donc un système d’équations linéaires lorsque l’on applique l’équation (42) au centre de chaque panneau, c’est à dire lorsque le point P devient le point de contrôle i avec i = 1,…,N. La N+1^è équation nécessaire provient de la condition de Kutta au bord de fuite. Le système d’équations ainsi obtenu peut prendre la forme matricielle suivante :

où [A] est la matrice des coefficients d’influence. Elle est obtenue à l’aide des équations précédentes et des conditions suivantes :

On résout le système de l’équation (49) en utilisant l’équation V_t,i = g_ipour obtenir directement les vitesses tangentielles. Puis on peut en déduire les coefficients de portance et de pression.

On définit que la vitesse tangentielle en un point du profil est égale à l’intensité du tourbillon. On va donc calculer la vitesse tangentielle au centre de chaque panneau, c’est-à-dire aux points de contrôles. Cependant, il est nécessaire de compléter ceci par la condition de Kutta. Aux bords de fuite:

tout simplement parce qu’on utilise la formulation par lignes de courant, sachant que sur le profil, on obtient :

avec L = r.V_¥.G selon le théorème de Kutta Joukowski, et q_¥ = (1/2).r.V_¥²

La circulation G peut s’écrire en fonction de l’intensité tourbillonnaire g de la manière suivante :

On peut maintenant approximer l’intégrale de contour par la sommation des longueurs des panneaux, ce qui amène à :

On obtient donc le coefficient de portance du profil à l‘aide de la formule suivante :

Remarque : On ne calcule pas le coefficient de traînée, car étant donné que le fluide étudié est considéré parfait donc non visqueux, la traînée occasionnée n’est pas représentative de la réalité.

Fonction menu ( )
Paramètres d’entrée : aucun
Paramètres de sortie : NACA, vitesse de l’écoulement, angle d’incidence, nombre de panneaux, entrés par l’utilisateur
Cette fonction va vous permettre de rentrer les paramètres de votre étude, comme le numéro de NACA, ou l’incidence du profil...

Fonction caracteristiques ( )
Paramètres d’entrée : NACA
Paramètres de sortie : symétrie, position de la cambrure maxi relative, cambrure maxi relative, épaisseur relative.
Cette fonction analyse le numéro de NACA pour en sortir les caractéristiques géométriques du profil.

Fonction discretisationNACA ( )
Paramètres d’entrée : symétrie, position de la cambrure maxi relative, cambrure maxi relative, épaisseur relative
Paramètres de sortie : vecteur des abscisses et vecteur des ordonnées des nœuds.
Cette fonction facétise le profil en panneaux de taille variable, et ajuste la densité de panneaux pour qu’il y en ait plus au niveau des bords d’attaque et de fuite (distribution cosinusoidale).

Fonction initmatrice ( )
Paramètres d’entrée : vecteur des abscisses et vecteur des ordonnées des nœuds.
Paramètres de sortie : vecteur des longueurs des panneaux, vecteur des abscisses et vecteur des ordonnées des points de contrôle, matrice A des coefficients d’influence.

Cette fonction permet de remplir la matrice du système linéaire de ses éléments.

Fonction init_b ( )
Paramètres d’entrée : vitesse de l’écoulement à l’infini, incidence, coordonnées de points de contrôles.
Paramètres de sortie : second membre initialisé

Fonction gauss_pivmax ( )
Paramètres d’entrée : la matrice du système, le vecteur second membre .
Paramètres de sortie : le vecteur des inconnues a l’emplacement mémoire du second membre (il l’écrase)

Cette fonction est un algorithme numérique de résolution des systèmes linéaires. Elle va nous permettre d’obtenir l’intensité du tourbillon sur chaque panneau.

7) Calcul des coefficients de pression et la composante du coefficient de portance sur chaque panneau du profil.

Fonction C_pression ( )
Paramètres d’entrée : le vecteur second membre, la vitesse à l’infini, le nombre de panneaux
Paramètres de sortie : le vecteur contenant tous les Cp aux panneaux

Paramètres d’entrée : le vecteur second membre, la vitesse à l’infini, le nombre de panneaux, la longueur des panneaux.
Paramètres de sortie : le vecteur des composantes de Cl a chaque panneau, et Le Cl du profil dans la dernière cas du vecteur.

Le but de ces fonctions est de calculer les coefficients de pression sur les panneaux et le coefficient de portance du profil.

Paramètres d’entrée : les paramètres a comparer
Paramètres de sortie : un fichier texte de résultats

Cette fonction nous permet de sortir les coordonnées des nœuds, points de contrôles, coefficients de pression, noms des auteurs, récapitulatif des différents paramètres du profil, coefficient de portance, afin d’être aisément lisibles pour un intervenant extérieur, et pourquoi pas, grâce à une extension du programme (non faite dans notre cas), permettre de venir prendre ces éléments pour refaire des calculs. Cette fonction sort de même un fichier texte nu (sans commentaires) avec les différents résultats, de façon a pouvoir les visualiser sous Excel ou tout autre traceur de graphes. Les noms des fichiers textes sont explicites. Resultataero.txt regroupe tous les résultats, sauf Cl en fonction de alpha.

Paramètres d’entrée : NACA
Paramètres de sortie : les nouveaux nombre de panneau, angle d’incidence de départ et d’arrivée, entrés par l’utilisateur.
Cette fonction va vous permettre de rentrer les paramètres des traitements supplémentaires, comme la nouvelle incidence du profil, le nouveau nombre de panneaux...

Paramètres d’entrée : le nombre d’éléments
Paramètres de sortie : la zone mémoire de la matrice
Cette fonction alloue une place mémoire en dynamique a une matrice, de façon à pouvoir libérer la mémoire grâce à free (nom de la matrice), ce qui optimise le code, et permet l’utilisation d’un grand nombre de panneaux (nous sommes arrivés jusqu’à 7000). Le calcul en statique nous aurait limité à 2 ou 300. dans notre cas, c’eut été suffisant, mais pas dans le cas d’une étude en liaison avec la ligne portante de prandtl.

Paramètres d’entrée : le nombre d’éléments
Paramètres de sortie : la zone mémoire du vecteur
La raison d’être de cette fonction est la même que précédemment, mais appliquée aux vecteurs, c’est à dire aux tableaux unidimensionnels.

Le programme a été réalisé en utilisant quelques astuces, comme par exemple stocker le coefficient de portance du profil à la fin du tableau des composantes de ce coefficient sur les panneaux du profil. Mais tout cela apparaît sur le listing clairement. L’algorithme numérique informel ne donnera pas accès à ces astuces. Il sera donc nécessaire de se reporter au programme joint au rapport.

Afin de vérifier l’exactitude des tracés de profils que nous avons obtenu, nous avons comparé nos résultats avec les tableaux de coordonnées de discrétisation fournis dans le livre de Ion Paraschivoiu. Comme on peut l’observer aux figures 5 et 6, les profils obtenus sont similaires aux profils théoriques.

Les résultats que nous obtenons pour les coefficients de pression sont parfaitement cohérents avec ceux fournis dans le livre de Ion Paraschivoiu.

Figure 7 Représentation du Cp en fonction de x/c pour un profil NACA 0018, d’angle d’incidence 5°

Figure 8 Représentation du Cp en fonction de l’angle d’incidence pour le profil NACA 0018

Ce graphique représente la variation du coefficient de pression sur l’intrados et l’extrados d’un profil symétrique NACA 0018 en fonction de la position sur le profil, en faisant varier l’incidence alpha du profil d’aile.

La courbe rouge au centre représente le profil sous incidence nulle : les pressions sont réparties de la même manière sur l’intrados et l’extrados, ce qui est normal. La vitesse tangentielle sur les parois du profil sont supérieures à la vitesse infini (sur le profil hors extrémités), car les lignes de courant du fluide, considéré comme parfait se rapprochent les unes des autres, ce qui, par application du théorème de conservation du débit, nous permet de dire que la vitesse le long du profil est supérieure. D’ou, par application du théorème de Bernoulli, on en déduit que la pression le long de ce même profil est inférieure à celle à l’infini. Ce phénomène est particulièrement visible juste après le bord d’attaque. Mais ce n’est pas vérifié au niveau du bord d’attaque lui même, qui est un point d’arrêt, donc où la vitesse s’ « annule », et ou le Cp tend vers 1. On a de même une décélération au bord de fuite qui fait augmenter le Cp.

L’augmentation de l’incidence permet de différencier les pressions à l’intrados et l’extrados. Sur l’extrados, la vitesse tangentielle augmente par rapport à la vitesse normale sous incidence nulle, donc, la pression diminue par rapport au profil sous incidence nulle.

À l’inverse, l’intrados voit sa vitesse tangentielle diminuer par rapport à précédemment, du fait de l’augmentation de l’obstacle que produit l’aile à l’écoulement venant de l’infini. D’où, augmentation de la pression par rapport au profil sous incidence nulle.

Le point d’arrêt est toujours le bord d’attaque, et la pression tend donc toujours vers 1 à ce niveau. Les courbes pour alpha 2, 4 et 5 en sont représentatives.

On constate à l’aide de la figure 9 que la variation du nombre de panneaux utilisés pour la discrétisation du profil, joue énormément sur le résultat du coefficient de portance

Figure 9 Influence du nombre de panneaux sur le coefficient de portance pour le NACA 0018

Figure 10 Représentation de l’évolution de la portance en fonction de l’angle d’incidence pour le profil NACA 0018

Ici, nous avons fait varier sur le programme l’angle d’incidence du profil symétrique NACA 0018, afin d’obtenir la portance du profil à diverses incidences. Nous remarquerons que Cl = f(alpha) suit une loi linéaire de pente environ égale à 0.12111. L’équation de Cl est Cl = 0.121*alpha pour le cas des profils symétriques 0018. C’est une relation simple, et donc aisée a manipuler. Ayant comparé pour plusieurs profils symétriques, nous nous rendons compte que l’équation de Cl sera alors toujours du type Cl = a*alpha où a est constant pour le profil. Cette relation nous permettra ainsi de connaître la portance quelle que soit l’incidence du profil symétrique, grâce, soit à l’essai sous une incidence quelconque, soit grâce à des tables qui donneraient directement l’équation de la droite.

Remarque : nous ne sommes pas allés au delà de 9 degrés, puisque la théorie fait fi de la viscosité, et de la couche limite, qui, au delà de 10 degrés d‘incidence, ne peut plus du tout être négligée : nos résultats ne seraient plus significatifs.

Figure 11 Influence du nombre de panneaux sur le résultat obtenu de la variation du coefficient de portance en fonction de l’incidence

Sur ce graphique, nous avons fait varier le nombre de panneaux de la discrétisation de 10 à 60. Le profil d’étude est encore symétrique de type NACA 0018. Pour ce faire, en nous basant sur l’étude précédente qui nous montre que la courbe obtenue est une droite passant par 0, nous avons uniquement cherché la valeur de la portance sous incidence de 9 degrés pour les différents nombres de panneaux.

Cette étude nous permet de montrer la convergence des résultats avec l’augmentation du nombre de panneaux sur le profil. La méthode des panneaux, de même que n’importe quelle méthode discrétionnaire, n’apporte pas une solution analytique au problème, mais une solution numérique. On pourrait croire que l’augmentation vers l’infini du nombre de panneaux permettrait d’obtenir un excellent résultat. C’est vrai sur ce point, mais c’est inutile, car ainsi que le montre la figure 11, nous avons une convergence rapide du résultat au delà de 40 panneaux.

Figure 12 Représentation de la variation du coefficient de portance en fonction de l’angle d’incidence pour le profil NACA 4412

La courbe précédente nous montre que pour les profils cambrés aussi, l’évolution du coefficient de portance du profil est une fonction linéaire de alpha, l’incidence du profil. Mais contrairement aux profils symétriques, les profils cambrés ont une portance initiale, ici de 0.37, ce qui fait que l’équation qui régit Cl(a) est : Cl(a) = a* a + b .

Cette relation nous permettra ainsi de connaître la portance quelle que soit l’incidence d’un profil cambré, grâce, soit à l’essai sous 2 incidences quelconques, soit grâce à des tables qui donneraient directement l’équation de la droite.

Les cinq méthodes de panneaux donnent de bons résultats, si on garde en mémoire qu’elles sont appliquées à des fluides parfaits, pour des incidences faibles. Parmi les méthodes les plus simples à calculer, il y a celle que nous avons étudiée : répartition de tourbillons par formulation de ligne de courant, ou encore celle du doublet linéaire.

Par contre, du fait même de ces limitations, il sera impossible d’étudier des profils dans des écoulement à forte viscosité, turbulents et sous forte incidence, ce qui nous laisse de nombreuses zones d’ombres dans l’étude d’une aile.

Néanmoins, il est possible d’apporter des coefficients correcteurs aux résultats, par exemple pour tenir compte de la couche limite, dans certaines méthodes, ce qui rapprochera encore les résultats numériques de la réalité.

Certaines méthodes de panneaux ne permettent pas non plus la simulation tridimensionnelle. C’est le cas de notre méthode par tourbillons et ligne de courant. De toutes les méthodes, celle qui a été l’objet de notre étude n’est donc pas à utilisée dans l’étude complète d’une aile.

On peut dire pour finir que ces méthodes ne nécessitent pas l’utilisation d’un nombre important de panneaux, et donc pas la résolution d’un système linéaire de très grande taille. Les calculs sont donc peu importants, ce qui rend cette méthode intéressante dans un but de design, afin de gagner du temps de conception.

De ce point de vue, notre méthode, ainsi que celle du doublet sont encore de loin les meilleures (à comparer avec Hess et Smith) .

PARASCHIVOIU, ION, Aérodynamique, Montréal, 2è édition, Ecole Polytechnique de Montréal, 1998, 572 pages.

IX. Listing du programme

visualisez le programme en C (bientôt le fichier .c ou .cpp en téléchargement, ainsi que le .exe)

Sommaire

1) Equations de base

2) Fondements théoriques

x_e = x – y_t.sinq

x_i = x + y_t.sinq

1) Méthode de Hess et Smith

2) Méthode avec tourbillons par formulation de vitesse

3) Méthode avec tourbillons par formulation de ligne de courant

4) Méthode avec doublets constants par formulation du potentiel

5) Méthode avec doublets linéaires par formulation du potentiel

1) Discrétisation

2) Equations de base

A_i,N+1 = -1

A_N+1,j = 0 pour j = 2,…,N-1

On résout le système de l’équation (49) en utilisant l’équation V_t,i = g_ipour obtenir directement les vitesses tangentielles. Puis on peut en déduire les coefficients de portance et de pression.

3) Conditions aux frontières:

4) Calcul de la distribution du gradient de pression.

5) Calcul du coefficient de portance.

1) Vérification de la discrétisation des profils

2) Résultats sur le coefficient de pression Cp

3) Résultats sur le coefficient de portance Cl

Sommaire

xe = x – yt.sinq

xi = x + yt.sinq

5) Méthode avec doublets linéaires par formulation du potentiel

Ai,N+1 = -1

AN+1,j = 0 pour j = 2,…,N-1

On résout le système de l’équation (49) en utilisant l’équation Vt,i = gi pour obtenir directement les vitesses tangentielles. Puis on peut en déduire les coefficients de portance et de pression.

x_e = x – y_t.sinq

x_i = x + y_t.sinq

A_i,N+1 = -1

A_N+1,j = 0 pour j = 2,…,N-1

On résout le système de l’équation (49) en utilisant l’équation V_t,i = g_ipour obtenir directement les vitesses tangentielles. Puis on peut en déduire les coefficients de portance et de pression.