FR 2952454 A1 20110513 FR 0957878 A 20091106 PROCEDE ET OUTIL DE SIMULATION DU COMPORTEMENT AERODYNAMIQUE D'UN ELEMENT AERODYNAMIQUE D'UN AERONEF PRESENTANT UN ANGLE DE FLECHE VARIABLE DESCRIPTION DOMAINE TECHNIQUE 10 La présente invention se rapporte au domaine général de l'aérodynamique et concerne la simulation numérique du comportement aérodynamique d'un élément aérodynamique d'un aéronef présentant un angle de flèche variable. 15 Elle trouve une application dans le domaine de l'aéronautique dans lequel la conception d'un aéronef nécessite de connaître précisément les coefficients aérodynamiques associés, par exemple, à son type de voilure. 20 ETAT DE LA TECHNIQUE ANTERIEURE Lors de la conception d'un aéronef et en particulier de la voilure de celui-ci, on cherche à déterminer les coefficients aérodynamiques globaux 25 associés à ladite voilure, tels que la portance, la traînée et le moment de tangage. Ces coefficients peuvent être déterminés de différentes manières, notamment par l'analyse d'un écoulement incident autour d'un profil d'aile 30 bidimensionnel de ladite voilure. Le profil d'aile est5 généralement obtenu à partir d'une coupe suivant un plan de l'une des ailes de ladite voilure, ledit plan du profil étant parallèle au plan de symétrie de l'aéronef. L'aéronef peut comporter une voilure présentant un angle de flèche cp constant. La figure 1 illustre, de manière schématique et en vue de dessus, une aile d'aéronef à angle de flèche constant dans le référentiel lié à l'aéronef. La figure 2 est une vue en perspective de l'aile représentée sur la figure 1. On définit un premier repère orthonormé direct Ra(Xa,Ya,Za), dit repère de l'aéronef, dont le plan (Xa,Za) coïncide avec le plan de symétrie de l'aéronef, et dont l'axe Ya est orthogonal au plan de symétrie de l'aéronef. Par souci de clarté et sans que cela soit limitatif, on considère ici que l'aile s'étend dans le plan (Xa, Ya) . Chacune des génératrices de l'aile forme un angle sensiblement constant par rapport à l'axe Ya. Cet angle est appelé angle de flèche. Ainsi, la génératrice de l'aile formant le bord d'attaque présente un angle de flèche de bord d'attaque (PBA. L'angle de bord de fuite cpBF est défini de la même manière, à partir de la génératrice de l'aile formant le bord de fuite. L'angle de flèche est dit constant lorsque les angles de flèche de bord d'attaque (PBA et de bord de fuite (pBF sont égaux. En d'autres termes, les génératrices de l'aile sont parallèles les unes aux autres. On peut définir un second repère orthonormé direct Rv(Xv,Yv,Zv), appelé repère local de la voilure, associé à chacune des génératrices de l'aile. Le plan (Xv,Yv) coïncide avec le plan (Xa, Ya) , et l'axe Yv est parallèle à l'axe de la génératrice considérée de la voilure. Ainsi, l'axe Xv du repère local de la voilure forme un angle local avec ledit plan (Xa,Za) égal à l'angle de flèche cp. Il est à noter que l'on passe du repère de l'aéronef Ra au repère local de la voilure Rv par une rotation élémentaire d'angle -cf) autour de l'axe Za. Il est à noter que la définition de l'angle de flèche n'est pas limitée à la voilure, mais peut également s'appliquer à d'autres éléments aérodynamiques de l'aéronef susceptibles de présenter un angle de flèche, comme, par exemple, les empennages. L'angle de flèche d'un empennage est alors défini à partir des génératrices de celui-ci, selon le même principe que présenté précédemment. Dans le référentiel lié à l'aéronef, l'aéronef est situé dans un écoulement tridimensionnel incident de vitesse amont, notée V00. Ladite vitesse amont V00 peut alors être décomposée, dans le repère local de la voilure Rv, en une première composante V0.coscp orientée dans le plan (Xv,Zv) et une seconde composante V00.sincp orientée dans le plan (Yv,Zv). Il est connu que la première composante V0.coscp détermine le champ des survitesses et donc la portance, alors que la seconde composante V00.sincp n'engendre aucune survitesse (P. Rebuffet, Aérodynamique expérimentale, 1950, p. 432-433). Pour déterminer les coefficients aérodynamiques, on néglige habituellement la seconde composante pour ne retenir que la première composante orientée dans le plan (Xv, Zv) . On considère alors un écoulement incident bidimensionnel autour d'un profil d'aile défini dans le plan du profil P, l'écoulement présentant une vitesse amont bidimensionnelle d'intensité V00.cos(p. Cependant, pour tenir compte du fait que le profil d'aile est défini par un plan du profil P parallèle au plan (Xa,Za) alors que la composante V0.cos(p de la vitesse amont tridimensionnelle V- est orientée dans le plan (Xv,Zv), le profil d'aile considéré est modifié. En effet, la composante V00.cos(p de la vitesse amont tridimensionnelle V- « voit », du fait de son orientation, un profil d'aile d'épaisseur relative plus importante que celle du profil défini dans le plan du profil P. En effet, l'épaisseur relative emax/c du profil défini dans le plan (Xv,Zv) est reliée à l'épaisseur relative du profil défini dans le plan du profil P par la relation suivante : 1 emax/ (1) x,~,Z) cos cl) C Âp) où emaX et c sont respectivement l'épaisseur maximale et la corde du profil dans le plan considéré, et (p est l'angle de flèche constant associé à ladite voilure. La corde est définie comme étant la distance, dans le plan considéré, entre le bord d'attaque et le bord de fuite. Aussi, le profil d'aile modifié est obtenu en corrigeant l'épaisseur locale du profil d'aile initial par la relation : e2 (x') = 1 .e1 (x') (2) cos (p où e1(x') et e2(x') sont les épaisseurs locales, respectivement, du profil initial et du profil modifié, à la distance x' au bord d'attaque suivant ladite corde définie dans le plan du profil P. La figure 3 illustre, en traits pointillés, le profil d'aile initial suivant le plan du profil P parallèle à (Xa,Za), et en trait plein, le profil d'aile modifié. L'étude de l'écoulement incident de vitesse amont bidimensionnelle V00.cos(p autour du profil d'aile modifié est habituellement réalisée par des essais en soufflerie ou par simulation numérique. A partir des champs de vitesse et de pression autour du profil obtenus par mesures expérimentales ou par la solution numérique, on peut calculer les coefficients locaux, tels que les coefficients de pression et de frottement pariétal, pour en déduire les coefficients aérodynamiques globaux mentionnés précédemment. Le comportement aérodynamique de l'aéronef est alors déterminé par lesdits coefficients aérodynamiques. Il est à noter que, dans le cas des essais en soufflerie, les coefficients aérodynamiques globaux peuvent être mesurés directement à l'aide de dispositifs de mesure d'efforts, tels que, par exemple, des balances trois ou six axes. Cependant, l'approche décrite précédemment ne peut simplement être étendue au cas des aéronefs dont la 30 voilure présente un angle de flèche variable. La figure 4 illustre, de manière schématique et en vue de dessus, une partie d'aile présentant un angle de flèche variable. L'angle de flèche de bord d'attaque (PBA est alors différent de l'angle de flèche de bord de fuite cpBF. De plus, pour chacune des génératrices de l'aile comprise entre le bord d'attaque et le bord de fuite, un angle de flèche local est associé qui présente une valeur différente. Ainsi, suivant la corde définie dans un plan du profil quelconque P parallèle au plan (Xa,Za), l'angle de flèche est défini localement et sa valeur varie entre une valeur de bord d'attaque et une valeur de bord de fuite selon une loi de variation donnée. A titre illustratif, dans le cas d'une aile trapézoïdale comme celle représentée sur la figure 4, la loi de variation de l'angle de flèche local est : tan «x') = tan (PBA. (1ù x') + tan (pBF .x' (3) Il est à noter que le repère local de la voilure Rv dépend ici de la valeur de l'angle de flèche local (p(x'). Les études en soufflerie et par simulation numérique réalisées dans le cas d'une voilure à angle de flèche constant sont difficiles à étendre au cas de la voilure à angle de flèche variable, dans la mesure où l'écoulement bidimensionnel dans cette dernière situation est accéléré, voire ralenti. En effet, l'angle de flèche de bord d'attaque (PBA est habituellement supérieur à l'angle de flèche de bord de fuite cpBF. Aussi, l'écoulement bidimensionnel incident de vitesse amont bidimensionnelle V0.cos(p autour du profil d'aile présente une vitesse au bord d'attaque dont l'intensité VOO.cos(PBA est inférieure à celle Vc.cospBF de la vitesse au bord de fuite. L'écoulement bidimensionnel est bien accéléré le long du profil d'aile. Dans le cas contraire, l'écoulement est ralenti. Or, les études bidimensionnelles réalisées en soufflerie ou par simulations numériques ne permettent actuellement pas de reproduire cette accélération ou ce ralentissement, ce qui nuit à la précision des résultats. Bien entendu, cette caractéristique de l'écoulement n'est présente que pour l'écoulement bidimensionnel de vitesse V0.cosp. L'écoulement tridimensionnel autour de l'aile tridimensionnelle ne presente pas cette caractéristique. On peut donc réaliser des études portant sur, non plus un écoulement simplifié bidimensionnel autour d'un profil d'aile, mais un écoulement tridimensionnel de vitesse amont tridimensionnelle V- autour d'une aile tridimensionnelle telle que représentée sur la figure 4. Ces études tridimensionnelles peuvent être réalisées en soufflerie ou par simulation numérique. Cependant, de part leur aspect tridimensionnel, la mise en oeuvre de ces études, tant sur le plan expérimental que numérique est rendue particulièrement lourde et complexe. Pour la même raison, l'analyse des données obtenues lors de ces études (par exemple, les champs de vitesse et de pression) est particulièrement longue et difficile. EXPOSÉ DE L'INVENTION L'invention a principalement pour but de présenter un procédé de simulation par ordinateur permettant d'obtenir le comportement aérodynamique d'un élément aérodynamique d'un aéronef, ledit élément aérodynamique présentant, entre le bord d'attaque et le bord de fuite de celui-ci, un angle de flèche local dont la valeur varie entre une valeur de bord d'attaque et une valeur de bord de fuite. Selon l'invention, ledit procédé comprend les étapes selon lesquelles : on obtient un profil bidimensionnel à partir d'une coupe suivant un plan d'un modèle géométrique tridimensionnel dudit élément aérodynamique ; - on réalise un maillage d'un domaine géométrique bidimensionnel délimité au moins en partie par ledit profil, ledit maillage définissant un domaine de calcul ; - on associe, à chaque maille dudit maillage, un repère tridimensionnel local R(X,Y,Z) pour lequel l'axe X forme un angle local avec ledit plan du profil, et l'axe Z est compris dans ledit plan du profil et orthogonal à la corde du profil, de sorte que, pour toute maille située entre le bord d'attaque et le bord de fuite du profil, l'angle local est égal à ladite valeur de l'angle de flèche localement définie ; - on résout par ordinateur, sur ledit maillage, un modèle numérique discret des équations de Navier-Stokes, de manière à obtenir la solution numérique d'un écoulement de fluide à l'intérieur dudit domaine de calcul, 8 ladite résolution dudit modèle numérique comportant, à chaque pas de temps, une étape de calcul de la vitesse dans chacune des mailles dudit maillage pour obtenir un champ de vitesse dans le domaine de calcul, ladite vitesse dans une maille considérée, exprimée dans le plan (X,Z) dudit repère tridimensionnel local correspondant, étant calculée à partir de la vitesse d'au moins une maille voisine amont ou aval suivant l'axe de ladite corde, ladite vitesse de ladite maille voisine étant préalablement exprimée dans le plan (X,Z) du repère tridimensionnel local de la maille considérée. Ledit plan du profil peut être choisi parallèle au 15 plan de symétrie dudit aéronef. Avantageusement, pour toute maille située en amont du bord d'attaque du profil, ledit angle local est égal à ladite valeur de l'angle de flèche au bord d'attaque ((PBA), et, pour toute maille située en aval du bord de 20 fuite du profil, ledit angle local est égal à ladite valeur de l'angle de flèche au bord de fuite ((PBF). De préférence, ledit modèle numérique discret des équations de Navier-Stokes comprend une condition imposée aux limites dudit maillage comportant une 25 vitesse locale définie, en chacune des mailles desdites limites, par la projection de la vitesse amont incidente dudit écoulement (V00) dans le plan (X,Z) dudit repère tridimensionnel local de la maille considérée. 30 Avantageusement, l'étape de calcul de la vitesse dans chacune des mailles dudit maillage comprend les sous-étapes suivantes . on définit une vitesse tridimensionnelle dans ladite maille voisine à la maille considérée à partir, d'une part de la vitesse bidimensionnelle dans ladite maille voisine, celle-ci étant exprimée dans le plan (X,Z) dudit repère tridimensionnel local correspondant, et d'autre part d'une troisième composante définie par la projection de la vitesse amont incidente dudit écoulement suivant l'axe Y du repère tridimensionnel local correspondant ; - on exprime ladite vitesse tridimensionnelle ainsi définie dans le repère tridimensionnel local de la maille considérée ; on projette ladite vitesse tridimensionnelle ainsi exprimée dans le plan (X,Z) du repère tridimensionnel local de la maille considérée, pour obtenir ainsi la vitesse bidimensionnelle dans ladite maille voisine exprimée dans le plan (X,Z) du repère tridimensionnel local de la maille considérée. De préférence, l'étape d'obtention dudit profil comprend une modification de l'épaisseur du profil suivant la relation : e2 (x') = 1 •el(x') cos(p(x ) où e1(x') et e2(x') sont les épaisseurs locales, respectivement, du profil initial et du profil modifié, à la distance x' au bord d'attaque suivant ladite corde, et cp(x') est la valeur de l'angle de flèche local à ladite distance x'. Avantageusement, ladite condition aux limites de vitesse imposée comporte, en outre, une composante locale de vitesse issue d'une perturbation du fluide environnant l'aéronef induite par ledit aéronef. Ladite solution numérique obtenue peut comprendre des champs de vitesse et de pression définis à l'intérieur du domaine de calcul. De préférence, une étape supplémentaire d'affichage de la solution numérique obtenue est prévue. L'invention porte également sur un procédé d'estimation de coefficients aérodynamiques d'un élément aérodynamique d'un aéronef dans un écoulement, ledit élément aérodynamique présentant, entre le bord d'attaque et le bord de fuite de celui-ci, un angle de flèche local dont la valeur varie entre une valeur de bord d'attaque et une valeur de bord de fuite. Selon l'invention, ce procédé comprend les étapes de . mise en oeuvre du procédé de simulation par ordinateur du comportement aérodynamique dudit élément aérodynamique selon l'une quelconque des caractéristiques précédentes ; puis - estimation des coefficients aérodynamiques dudit aéronef à partir de ladite solution numérique obtenue. Ladite étape d'estimation des coefficients aérodynamiques peut comprendre une étape intermédiaire d'estimation des coefficients locaux de pression et de frottement pariétal à partir de ladite solution numérique obtenue et de la pression dynamique définie à partir de la vitesse amont incidente dudit écoulement. Lesdits coefficients aérodynamiques locaux peuvent être intégrés suivant ledit profil initial à l'épaisseur non modifiée, de manière à obtenir des coefficients aérodynamiques globaux exprimés dans le repère lié à l'aéronef. Lesdits coefficients aérodynamiques globaux ainsi obtenus peuvent être projetés dans un référentiel aérodynamique orienté suivant ladite vitesse moyenne amont, de manière à obtenir les coefficients de portance, de traînée et du moment de tangage. D'autres avantages et caractéristiques de l'invention apparaîtront dans la description détaillée non limitative ci-dessous. BRÈVE DESCRIPTION DES DESSINS On décrira à présent, à titre d'exemples non limitatifs, des modes de réalisation de l'invention, en se référant aux dessins annexés, dans lesquels : La figure 1, déjà décrite, est une vue schématique de dessus d'une aile d'aéronef présentant un angle de flèche constant ; La figure 2, déjà décrite, est une vue schématique en perspective de l'aile représentée sur la figure 1 ; La figure 3, déjà décrite, est une vue schématique de profils d'aile non modifiés et modifiés obtenus à partir de la voilure représentée sur les figures 1 et 2 ; La figure 4, déjà décrite, est une vue schématique de dessus d'une aile de voilure présentant un angle de flèche variable ; La figure 5 est une vue schématique du domaine 5 géométrique bidimensionnel correspondant au domaine de calcul ; La figure 6 est une représentation schématique de quelques mailles appartenant au maillage, dans le cas d'un maillage cartésien régulier ; et 10 La figure 7 est une vue schématique de dessus d'une partie d'aile sur laquelle sont représentés quelques repères tridimensionnels locaux associés à différentes mailles du maillage. 15 EXPOSÉ DÉTAILLÉ D'UN MODE DE RÉALISATION PREFERE Les modes de réalisation de l'invention qui sont présentés par la suite portent sur une voilure d'un aéronef dont les ailes présentent un angle de flèche variable. Cependant, la présente invention n'est pas 20 limitée aux ailes de la voilure de l'aéronef. Tout élément aérodynamique de l'aéronef peut être étudié, comme, par exemple, les empennages. On considère un aéronef dont la voilure, représentée sur la figure 4, présente un angle de 25 flèche variable. Plus précisément, ladite voilure comporte deux ailes dont chacune présente, entre son bord d'attaque et son bord de fuite, un angle de flèche local dont la valeur varie entre une valeur de bord d'attaque et une valeur de bord de fuite selon une loi 30 de variation donnée. Le bord d'attaque et le bord de fuite sont définis par la génératrice correspondante de l'aile. La corde est la distance reliant le bord d'attaque et le bord de fuite suivant un plan de coupe transversale de l'aile, de préférence parallèle au plan de symétrie de l'aéronef. Dans le référentiel lié à l'aéronef, un écoulement incident tridimensionnel, de vecteur vitesse amont V00, s'écoule autour de l'aéronef. On considère un repère de l'aéronef Ra et un repère local de la voilure Rv tels que définis précédemment. On définit également un repère orthonormé direct Re(Xe,Ye,Ze), dit repère aérodynamique de l'écoulement incident, dont l'axe Xe des abscisses est orienté suivant la direction du vecteur vitesse amont V- de l'écoulement incident et dont l'axe Ye est orthogonal au plan de symétrie de l'aéronef. Le vecteur vitesse amont V- peut former un angle non nul, appelé ici angle d'incidence î, avec l'axe Xa du repère de l'aéronef. Cet angle est parfois appelé angle d'attaque (AoA, en anglais). Par la suite, et dans un souci de clarté de la description sans que cela ne soit limitatif, on considère que le vecteur vitesse amont V- de l'écoulement incident est parallèle au plan de symétrie de l'aéronef. En d'autres termes, l'aile présente un angle de dérapage sensiblement nul. Les plans (Xe,Ze) et (Xa,Za) sont alors parallèles. Il est à noter que le passage du repère aérodynamique au repère de l'aéronef s'effectue par rotation autour de l'axe Ya d'un angle égal à l'angle d'incidence -a. Le repère aérodynamique et le repère de l'aéronef coïncident donc l'un avec l'autre lorsque l'angle d'incidence a est nul. Par ailleurs, le passage du repère de l'aéronef Ra au repère local de la voilure Rv s'effectue par rotation autour de l'axe Za d'un angle égal à l'angle de flèche local -p(x). Le procédé de simulation par ordinateur du comportement aérodynamique dudit élément aérodynamique, ici de la voilure de l'aéronef, selon l'invention est maintenant décrit, en référence aux figures 5 à 7. A partir d'un modèle géométrique tridimensionnel de la voilure dudit aéronef, on réalise un profil d'aile bidimensionnel 1. Le profil 1 est obtenu à partir d'une coupe suivant un plan, appelé plan du profil P, du modèle géométrique de la voilure. Le plan du profil P est parallèle au plan de symétrie de l'aéronef, donc au plan (Xa,Za) du repère de l'aéronef. Comme expliqué précédemment, l'épaisseur locale du profil d'aile 1 est modifiée. Cette modification provient du fait que l'on considère, comme explicité précédemment, le long du profil 1 défini dans le plan du profil P, un écoulement bidimensionnel dont la vitesse est la composante V00.cosp selon le plan (Xv,Zv) de la vitesse amont tridimensionnelle Vo0. Or, le plan (Xv,Zv) et le plan du profil P ne coïncident pas. De ce fait, la composante Vo.cosp de la vitesse amont tridimensionnelle Vo0 « voit » un profil différent du profil 1 considéré, à savoir le profil défini selon le plan (Xv,Zv). Aussi, l'épaisseur locale du profil d'aile considéré, à savoir le profil d'aile 1 défini selon le plan du profil P, est modifiée suivant la relation suivante : = 1 .e (x » cos (p (x') où e1(x') et e2(x') sont les épaisseurs locales, respectivement, du profil initial et du profil modifié, à la distance x' au bord d'attaque suivant la corde du profil, et cp(x') est la valeur de l'angle de flèche 10 local à ladite distance x'. Il est à noter que l'angle de flèche considéré pour réaliser cette modification de l'épaisseur locale est bien l'angle de flèche local et non pas une valeur de l'angle de flèche qui reste constante quelle que soit 15 la distance x' suivant la corde du profil. Un domaine géométrique bidimensionnel 2 est défini, qui correspond à la zone d'espace qui entoure le profil d'aile 1 et dans laquelle l'écoulement de l'air sera 20 reproduit par simulation. Ce domaine 2 est également appelé domaine de calcul. Ce domaine géométrique 2 présente, de préférence, une forme de rectangle, mais peut présenter toute autre forme, par exemple un rectangle dont une face est 25 convexe lorsqu'elle est vue de l'extérieur du rectangle. La face amont 3 du domaine géométrique est située face au nez du profil d'aile et la face aval 4 est disposée en vis-à-vis de la face amont. Les faces amont (4) 3 et aval 4 sont reliées par deux faces inférieure 5 et supérieure 6. Le profil d'aile 1 est disposé dans ce domaine géométrique 2 de manière à être suffisamment éloigné de chacune desdites faces. Les faces amont 3 et aval 4 peuvent être disposées à une distance du profil d'aile 1 de l'ordre d'une à plusieurs dizaines de fois la corde du profil. Les faces inférieure 5 et supérieure 6 peuvent être situées à une distance du profil d'aile 1 de l'ordre de une à plusieurs fois la corde du profil. Un modèle mathématique est choisi pour reproduire l'écoulement de l'air environnant le profil d'aile. Ce modèle est le modèle des fluides réels compressibles, dit de Navier-Stokes. Le fluide est supposé visqueux, newtonien et compressible, et conducteur de chaleur. Ce modèle physique, bien connu, comprend un ensemble d'équations aux dérivées partielles appelées équation de continuité, équation de la dynamique et équation de l'énergie et s'écrivant, respectivement : dp +p aU'=0 dt axe dU~ aP a2U. 1 a 'aUz p =più +µ + dt axe aU axe axe 3 axe axl aaT pC, =ùP '+(D, dt axe axe axe où p (x) la masse volumique du fluide, U (x) le champ de vitesse, P (x) le champ de pression, T (x) le champ de température, F est une force extérieure de volume, comme la gravité, µ la viscosité dynamique, CV la chaleur spécifique à volume constant, d)ä le taux de dissipation par viscosité, et a, la conductivité thermique. Ce modèle physique de Navier-Stokes peut être cependant adapté à la description des écoulements turbulents. Les équations générales de Navier-Stokes données ci-dessus peuvent être reformulées sous la forme connue des équations moyennées, appelées RANS pour Reynolds Averaged Navier-Stokes, et comprennent alors un modèle de turbulence. Différents modèles connus de turbulence peuvent être utilisés, comme, par exemple, les modèles de turbulence à deux équations tels que les modèles k-E et k-w. Parmi ces modèles, le modèle SST tel que décrit dans l'article de Menter intitulé « Zonal Two Equation k-w Turbulence Models for Aerodynamic Flows » et publié en 1993 dans AIAA Paper 93-2906 peut être utilisé. La modification dite de Kato-Launder peut être prise en compte pour corriger le terme de production turbulente dans l'équation de l'énergie cinétique turbulente. Cette approche est décrite dans l'article de Kato et Launder intitulé « The Modeling of Turbulent Flow Around Stationary and Vibrating Square Cylindres » et publié en 1993 dans Proc. 9th Symposium on Turbulent Shear Flows, Kyoto, pages 10.4.1 à 10.4.6. Par ailleurs, il est également possible de prendre en compte, dans le modèle SST avec la modification Kato-Launder, la correction dite Kok TNT qui diminue la surproduction turbulente au coeur des tourbillons. Ce modèle est décrit dans l'article de Kok intitulé « Resolving the dependence of free stream values of the k-w turbulence model » et publié en 2000 dans AIAA Journal 38, 1292-1295. Bien entendu, l'invention n'est pas limitée à l'utilisation de ces modèles de turbulence. D'autres modèles peuvent être utilisés, comme par exemple le modèle k-w EARSM développé par Hellsten. Le modèle physique de Navier-Stokes comprend également une condition initiale et des conditions aux limites. De manière à reproduire l'écoulement accéléré présenté précédemment, on étend, à chaque point du domaine géométrique bidimensionnel, le repère tridimensionnel local Rv(Xv,Yv,Zv) de la voilure en définissant un repère tridimensionnel local R(X,Y,Z), celui-ci étant de préférence orthogonal. L'axe X du repère tridimensionnel local R forme un angle local cp (x, z) avec le plan du profil P et l'axe Z est compris dans ledit plan du profil et est orthogonal à ladite corde. L'angle local cp(x,z) est défini en tout point de coordonnées (x,z) du domaine géométrique bidimensionnel. Cet angle local cp(x,z) peut être défini de la manière suivante : - Pour tout point situé en amont du bord d'attaque du profil suivant l'axe de la corde de celui-ci, l'angle local cp est égal à ladite valeur de l'angle de flèche cp au bord d'attaque : bx _xBF, Vz, (p (x, z) =(PBF, xBF est l'abscisse du bord de fuite. Ainsi, en chaque point du domaine géométrique bidimensionnel 2, un repère tridimensionnel local est défini qui dépend de l'angle de flèche local. Le domaine géométrique bidimensionnel 2 est ainsi formé de trois parties Al, A2 et A3, comme l'illustre la figure 5. Il est à noter que l'angle local cp(x,z) est constant dans les parties amont Al et aval A3. Dans la partie A2, il coïncide avec l'angle de flèche local et ne dépend que de la position en x. Ainsi, dans la partie A2, le repère tridimensionnel local R coïncide avec le repère local de la voilure Rv. La condition initiale peut comprendre un champ de vitesse imposé dans tout le domaine géométrique bidimensionnel. La vitesse à l'instant initial, en chacun des 30 points du domaine géométrique bidimensionnel, peut être la projection de la vitesse tridimensionnelle amont V-dans le plan (X,Z) du repère tridimensionnel local R correspondant. A titre illustratif, dans le cas d'une vitesse amont V- présentant un angle d'incidence î nul, la projection de la vitesse amont dans le plan (X,Z) du repère tridimensionnel local correspondant s'écrit : (Voo. cos [(p (x, z) ] , 0) , où l'angle local (p (x, z) dépend de la position du point considéré selon les trois parties géométrique bidimensionnel définies Ainsi, la condition initiale correspond un écoulement accéléré lorsque mesure où la vitesse bidimensionnelle = Voo. cos(pBF en aval du profil d'aile la vitesse bidimensionnelle V- . cos [cf) (x Les conditions aux limites des différentes frontières (faces du domaine et surface du profil 20 d'aile) du domaine géométrique peuvent se répartir en plusieurs catégories. Une première catégorie de conditions aux limites concerne la condition imposée à la surface du profil d'aile. 25 La condition aux limites imposée à la surface du profil d'aile est une condition de type « paroi » classique. Par condition de type « paroi » classique, on entend que la vitesse à la surface considérée est 30 nulle, plus particulièrement les composantes tangentielle et normale, et que toute perturbation du du domaine précédemment. effectivement à (PBA > cpBF, dans la V-. cos [(p (x>_xBF, z) ] est supérieure à fluide environnant l'aéronef est réfléchie par ladite surface. Une seconde catégorie de conditions aux limites concerne les faces amont 3, aval 4, inférieure 5 et supérieure 6 du domaine géométrique 2. Concernant la face amont 3, une condition de vitesse est imposée qui présente une composante fixe obtenue par la projection de la vitesse tridimensionnelle amont V- dans le plan (X,Z) du repère tridimensionnel local R au point considéré appartenant à ladite face amont 3. La vitesse imposée comprend avantageusement une composante supplémentaire V' issue d'une perturbation du fluide environnant l'aéronef induite par ledit aéronef. Ainsi, un flux d'air entre dans le domaine géométrique avec une vitesse moyenne correspondant à Vco.coscpBA lorsque l'angle d'incidence est sensiblement nul. Une perturbation du fluide environnant peut être tout type d'onde de pression telle qu'une onde sonore, une onde de choc ou plus généralement une onde de discontinuité, voire également un flux d'air induit par la présence du profil d'aile. La perturbation de type flux d'air peut provenir directement du profil d'aile. De plus, toute perturbation du fluide environnant induite par l'aéronef peut sortir du domaine géométrique 2 au travers de la face amont 3. La perturbation est alors dite sortante. La vitesse imposée à la face amont 3 peut être non uniforme et dépend de la perturbation sortante. Ainsi, la composante de vitesse Vco.cos(pBA est de préférence constante quel que soit le point considéré de la face amont 3, alors que la nécessairement constante, en orientation, et dépend du point amont 3. Il est à noter que ladite perturbation sortante peut être sensiblement nulle ou négligeable devant le flux d'air entrant au travers de la face amont. La vitesse imposée est alors sensiblement égale à la composante V00. cos(PBA. Concernant les autres faces du domaine géométrique, c'est-à-dire les faces aval 4, inférieure 5 et supérieure 6, la condition aux limites imposée peut être du même type que pour la face amont. A chacune des faces aval 4, inférieure 5 et supérieure 6, la condition aux limites présente une vitesse imposée qui présente une composante fixe obtenue par la projection de la vitesse tridimensionnelle amont dans le plan (X,Z) du repère tridimensionnel local R au point considéré appartenant à ladite face considérée. La vitesse imposée comprend avantageusement une composante supplémentaire V' issue d'une perturbation du fluide environnant l'aéronef induite par ledit aéronef. Toute perturbation du fluide environnant induite 25 par l'aéronef peut alors sortir du domaine géométrique au travers desdites faces. Comme précédemment, la composante supplémentaire de la vitesse imposée sur ces faces n'est pas nécessairement constante, en intensité et en 30 orientation, selon le point de la face considérée. vitesse V' n'est pas intensité et en considéré de la face A titre illustratif, dans le cas d'une vitesse amont tridimensionnelle V- présentant un angle d'incidence î nul, la vitesse imposée suivant la face amont est (Voo. cos(pBA, 0) +V' , suivant la face aval (Voo. cos(pBF, 0) +V' , et suivant les faces inférieure et supérieure : (V0. cos [(p (x) ] , 0) +V' . Ainsi, le modèle physique comprend, comme décrit précédemment, l'ensemble d'équations du modèle de Navier-Stokes, de préférence adapté à la description d'un écoulement turbulent, ainsi qu'une condition initiale, et des conditions aux limites imposées aux frontières du domaine géométrique (faces du domaine et surface du modèle d'aéronef). Le domaine géométrique est discrétisé pour obtenir un maillage bidimensionnel. Les frontières du maillage coïncident avec celles du domaine géométrique. Le maillage peut être du type structuré ou non structuré. Le maillage est réalisé à l'aide d'un logiciel, par exemple le logiciel CATIA V5. Dans le cas général, chaque maille est numérotée en i, k et correspond à une unique position définie en x, z dans le domaine géométrique bidimensionnel. Ainsi, les coordonnées x et z sont chacune fonctions de i et de k : x (i, k) et z (i, k) . La figure 6 illustre une partie d'un exemple de maillage cartésien non orthogonal à face centrée. Il est à noter que le vecteur i n'est pas nécessairement colinéaire à l'axe de la corde. Le maillage comprend les trois parties Al, A2 et A3 du domaine géométrique bidimensionnel. Aussi, le repère tridimensionnel local R dépend de la valeur de l'angle local (p (i, k) . Ainsi . - Pour toute maille située en amont du bord d'attaque du profil suivant l'axe de la corde de celui-ci, l'angle local (p est de préférence égal à ladite valeur de l'angle de flèche (p au bord d'attaque : b(i,k)/x(i,k) Un modèle numérique discret est ensuite obtenu par la discrétisation temporelle et spatiale du modèle physique de Navier-Stokes décrit précédemment. Différents types de schémas numériques généraux peuvent être utilisés pour la discrétisation spatiale, tels que les volumes finis, les éléments finis ou les différences finies pour la discrétisation spatiale. De préférence, la technique des volumes finis est utilisée, notamment du deuxième ordre à faces centrées, et la discrétisation en temps est obtenue par un schéma explicite de Runge-Kutta. La discrétisation temporelle peut également être obtenue par divers schémas connus de l'homme du métier. Ainsi, on obtient un modèle numérique discret de Navier-Stokes comprenant une condition initiale imposée au temps initial à toutes les mailles du domaine géométrique, et des conditions aux limites imposées aux frontières du maillage. La résolution du modèle numérique discret décrit précédemment est effectuée par ordinateur, ce qui permet de simuler l'écoulement de l'air autour du profil d'aile modifié dudit aéronef. Selon l'invention, ladite résolution du modèle numérique comporte, à chaque pas de temps, une étape de calcul de la vitesse dans chacune des mailles dudit maillage pour obtenir un champ de vitesse dans le domaine de calcul, ladite vitesse dans une maille considérée, exprimée dans le plan (X,Z) dudit repère tridimensionnel local correspondant, étant calculée à partir de la vitesse d'au moins une maille voisine amont ou aval suivant l'axe de ladite corde, ladite vitesse de ladite maille voisine étant préalablement exprimée dans le plan (X,Z) du repère tridimensionnel local de la maille considérée. Aux mailles (i, k) , (i-1, k) et (i+l, k) sont associés, respectivement, les repères tridimensionnels locaux R(i'», R(i-1'» et R(i"'» Comme l'illustre la figure 7, la vitesse V(i'k) de la maille (i,k) est exprimée dans le plan (X'i'') , Z (i' ) du repère tridimensionnel local R(i'k). Pour souci de clarté, l'indice k n'est pas indiqué sur cette figure 7. Les vitesses V(+1'k) de la maille voisine (i-1, k) suivant l'axe de la corde et la vitesse V(i+1'k) de la maille voisine (i+l,k) sont exprimées, respectivement, dans les plans (X (i-1' k) , Z (i-1' k)) et (X (++1, k) , z (++1, k)) correspondant aux repères tridimensionnels locaux R(+-1,k) et R(i+1,k) La vitesse V(i'k) de la maille (i,k) est calculée, par exemple, à partir des vitesses V(i-1,k) et V(i+1'k) Cependant, les vitesses V(i-1,k) et V(i+1,k) utilisées pour cette étape de calcul sont chacune préalablement projetées dans le repère tridimensionnel local R(i'k) associé à la maille considérée (i,k). Il est à noter que l'on passe du repère R(+"k), i'i, au repère R(i'k) par une simple rotation autour de l'axe Z d'un angle [''p' k) i Par ailleurs, il est à noter que, par maille voisine amont ou aval suivant l'axe de la corde, on entend toute maille autre que la maille considérée située en amont ou en aval de celle-ci suivant l'axe de la corde. La maille voisine peut être une maille jouxtant la maille considérée. Elle peut également être disposée à distance de ladite maille considérée. Ainsi, pour obtenir la vitesse dans une maille considérée, toutes les vitesses utilisées correspondant aux mailles voisines amont ou aval suivant l'axe de la corde sont projetées dans le repère tridimensionnel local R(i'k) associé à la maille considérée (i, k) . Ainsi, sans avoir recours à des sources ou des puits de fluide, l'écoulement est accéléré (dans le cas présent où (PBA > (pBF). Bien entendu, dans le cas où (PBA Ainsi, dans le cas où l'angle d'incidence est nul, la vitesse bidimensionnelle V2D~++1,10 dans ladite maille voisine (i+l,k) s'écrit (vxi+1, k, vzi+1, k) dans le plan (X(i+1,k),Zci+1,k)) dudit repère tridimensionnel local correspondant R(i+1,k) La troisième composante est définie par la projection de la vitesse amont tridimensionnelle suivant l'axe Y(i+1,k) du repère tridimensionnel local R(i+1,k) correspondant et s'écrit V-.sin[cp(i+l,k)]. La vitesse tridimensionnelle V3D(++1, k) dans ladite maille voisine s'écrit alors (Vxi+1,k, Vue. sin [(p (i+1, k) ] , Vzi+1,k) dans le repère tridimensionnel local correspondant R(i+1,k) Ensuite, on exprime ladite vitesse tridimensionnelle V3D(i+i,k) ainsi définie dans le repère tridimensionnel local R(+'k) de la maille considérée. On note V3D(i+1,k) IR(i+1,k) la vitesse tridimensionnelle définie dans le repère R(i+1,k), et V3D(i+i,k) IR(i,k) la même vitesse tridimensionnelle exprimée dans le repère R(i'k) On a alors la relation de passage entre les deux repères . (i+1, k) (i+1, k)R(i, k)] U3D(i+1, k) R(1,k) V3D R(1+1, k) = PLR où P [R(i+1,k)R(i'k) ] est la matrice de passage du repère R(i+1'k) au repère R(i'k) Enfin, on projette ladite vitesse tridimensionnelle V3D (i+1, k) I R(i, k) ainsi exprimée dans le plan (X', Z (i, k) ) du repère tridimensionnel local R(i'k) de la maille considérée, pour obtenir ainsi la vitesse bidimensionnelle V2D(i+1,k) IR(i,k) dans ladite maille voisine (i+l, k) exprimée dans le plan (X(i'k), Z(i'k)) du repère tridimensionnel local R(i'k) de la maille considérée. Cette étape de calcul du champ de vitesse est réalisée à chaque pas de temps, pour toutes les mailles du maillage. Ainsi, sans avoir recours à des sources ou des puits de fluide, l'écoulement est accéléré (dans le cas présent où (PBA > (pBF). Bien entendu, dans le cas où (PBA Enfin, pour suivre l'évolution de la résolution du modèle numérique discret, un critère de convergence est utilisé qui permet d'arrêter la simulation lorsqu'il est vérifié. Le critère de convergence peut être une grandeur physique telle qu'une vitesse ou une pression mesurée en un point donné du maillage. Lorsque cette donnée est stationnaire, on considère que l'écoulement du fluide dans le domaine de calcul est établi. La simulation peut alors être stoppée. On obtient une solution numérique de l'écoulement du fluide environnant le profil d'aile. Cette solution numérique comprend des champs de vitesse, de pression, de température en tout point du maillage. Elle est stockée dans la mémoire dudit ordinateur, ou dans des moyens de stockage séparés. La solution numérique peut être affichée à l'écran d'un ordinateur. Par des moyens logiciels, le comportement aérodynamique de l'aéronef est analysé à partir de la solution numérique obtenue. On calcule alors, à partir des champs de vitesse et de pression obtenus, les coefficients aérodynamiques locaux du profil d'aile modifié de l'aéronef. Les coefficients aérodynamiques locaux comprennent notamment les coefficients de pression et de frottement. Pour l'adimensionnement, on utilise, de préférence, comme pression de référence la pression dynamique tridimensionnelle, et non pas la pression dynamique bidimensionnelle. En effet, celle-ci n'est pas définie de manière unique puisqu'elle peut varier dans le domaine de calcul. Pour obtenir les coefficients aérodynamiques globaux dans le repère de l'aéronef, on intègre les coefficients locaux précédemment obtenus le long du profil d'aile. Pour l'adimensionnement des longueurs, on peut utiliser la corde du profil. On obtient alors les coefficients globaux dans le repère de l'aéronef, à savoir les coefficients de force normale, de force tangentielle et le coefficient de moment de tangage. Enfin, la projection dans un repère aérodynamique de l'écoulement Re des coefficients globaux précédemment calculés permet d'obtenir les coefficients aérodynamiques globaux dans ledit repère aérodynamique Re, c'est-à-dire les coefficients de portance, de traînée et de moment de tangage. La simulation peut être répétée pour différentes valeurs de la vitesse V00, de la loi de variation de l'angle de flèche p(x') suivant la corde du profil, ainsi que de l'angle a. Cette analyse peut également être effectuée pour divers modèles géométriques d'aéronef, et donc différents profils d'aile. Les résultats d'analyse peuvent être classés sous forme de bases de données. Chaque base de données indique les coefficients aérodynamiques de l'aéronef pour une large gamme du domaine de vol défini par, notamment, la vitesse de l'aéronef, la loi de variation de l'angle de flèche p(x'), l'angle î et le modèle d'aéronef utilisé. Bien entendu, diverses modifications peuvent être apportées par l'homme du métier à l'invention qui vient d'être décrite, uniquement à titre d'exemples non limitatifs. Concernant le modèle physique utilisé, on peut choisir un modèle simplifié de Navier-Stokes. Par exemple, le fluide environnant l'aéronef peut être considéré comme non visqueux. Le modèle de Navier-Stokes se réduit alors aux équations d'Euler pour les fluides parfaits. Par ailleurs, la description de la turbulence par le modèle de Navier-Stokes peut être obtenue, alternativement au modèle RANS décrit précédemment, par les modèles du type LES (Large Eddy Simulation), du type hybride RANS/LES, ou DES (Detached Eddy Simulation). La présente invention n'est pas limitée au cas de l'angle de dérapage nul. Dans le cas d'un angle de dérapage non nul, le passage du repère aérodynamique au repère de l'aéronef s'effectue par deux rotations élémentaires, une première rotation autour de l'axe Ya d'un angle égal à l'angle d'incidence oc et une seconde rotation autour de l'axe Za d'un angle égal à l'angle de dérapage P. Le procédé est ensuite similaire à celui décrit précédemment.20 REVENDICATIONS1. Procédé de simulation par ordinateur du comportement aérodynamique d'un élément aérodynamique d'un aéronef dans un écoulement, ledit élément aérodynamique présentant, entre le bord d'attaque et le bord de fuite de celui-ci, un angle de flèche local (q)) dont la valeur varie entre une valeur de bord d'attaque et une valeur de bord de fuite, caractérisé en ce que : - on obtient un profil bidimensionnel à partir d'une coupe suivant un plan (P) d'un modèle géométrique tridimensionnel dudit élément aérodynamique ; - on réalise un maillage d'un domaine géométrique bidimensionnel délimité au moins en partie par ledit profil, ledit maillage définissant un domaine de calcul ; - on associe, à chaque maille dudit maillage, un repère tridimensionnel local R(X,Y,Z) pour lequel l'axe X forme un angle local avec ledit plan du profil, et l'axe Z est compris dans ledit plan du profil et orthogonal à la corde du profil, de sorte que pour toute maille située entre le bord d'attaque et le bord de fuite du profil, l'angle local est égal à ladite valeur de l'angle de flèche (y) localement définie ; - on résout par ordinateur, sur ledit maillage, un modèle numérique discret des équations de Navier-Stokes, de manière à obtenir la solution numérique d'un écoulement de fluide à l'intérieur dudit domaine de calcul,ladite résolution dudit modèle numérique comportant, à chaque pas de temps, une étape de calcul de la vitesse dans chacune des mailles dudit maillage pour obtenir un champ de vitesse dans le domaine de calcul, ladite vitesse dans une maille considérée, exprimée dans le plan (X,Z) dudit repère tridimensionnel local correspondant, étant calculée à partir de la vitesse d'au moins une maille voisine amont ou aval suivant l'axe de ladite corde, ladite vitesse de ladite maille voisine étant préalablement exprimée dans le plan (X,Z) du repère tridimensionnel local de la maille considérée. 2. Procédé de simulation par ordinateur du comportement aérodynamique d'un élément aérodynamique d'un aéronef selon la revendication 1, caractérisé en ce que, pour toute maille située en amont du bord d'attaque du profil, ledit angle local est égal à ladite valeur de l'angle de flèche au bord d'attaque (qBA), et en ce que, pour toute maille située en aval du bord de fuite du profil, ledit angle local est égal à ladite valeur de l'angle de flèche au bord de fuite ((PBF) 3. Procédé de simulation par ordinateur du comportement aérodynamique d'un élément aérodynamique d'un aéronef selon la revendication 1 ou 2, caractérisé en ce que, ledit modèle numérique discret des équations de Navier-Stokes comprend une condition imposée aux limites dudit maillage comportant une vitesse locale définie, en chacune des mailles desdites limites, parla projection de la vitesse amont incidente dudit écoulement (Vao) dans le plan (X,Z) dudit repère tridimensionnel local de la maille considérée. 4. Procédé de simulation par ordinateur du comportement aérodynamique d'un élément aérodynamique d'un aéronef selon la revendication 3, caractérisé en ce que ladite condition aux limites de vitesse imposée comporte, en outre, une composante locale de vitesse issue d'une perturbation du fluide environnant l'aéronef induite par ledit aéronef. 5. Procédé de simulation par ordinateur du comportement aérodynamique d'un élément aérodynamique d'un aéronef selon l'une quelconque des revendications 1 à 4, caractérisé en ce que l'étape de calcul de la vitesse dans chacune des mailles dudit maillage comprend les sous-étapes suivantes : - on définit une vitesse tridimensionnelle dans ladite maille voisine à la maille considérée à partir, d'une part de la vitesse bidimensionnelle dans ladite maille voisine, celle-ci étant exprimée dans le plan (X,Z) dudit repère tridimensionnel local correspondant, et d'autre part d'une troisième composante définie par la projection de la vitesse amont incidente dudit écoulement suivant l'axe Y du repère tridimensionnel local correspondant ; - on exprime ladite vitesse tridimensionnelle ainsi définie dans le repère tridimensionnel local de la 30 maille considérée ;on projette ladite vitesse tridimensionnelle ainsi exprimée dans le plan (X,Z) du repère tridimensionnel local de la maille considérée, pour obtenir ainsi la vitesse bidimensionnelle dans ladite maille voisine exprimée dans le plan (X,Z) du repère tridimensionnel local de la maille considérée_ 6. Procédé de simulation par ordinateur du comportement aérodynamique d'un élément aérodynamique d'un aéronef selon l'une quelconque des revendications 1 à 5, caractérisé en ce que l'étape d'obtention dudit profil comprend une modification de l'épaisseur du profil suivant la relation : e2 (x') = 1 •e1(x') cos çp(x ) où e1(x') et e2(x') sont les épaisseurs locales, respectivement, du profil initial et du profil modifié, à la distance x' au bord d'attaque suivant ladite corde, et cp(x') est la valeur de l'angle de flèche local à ladite distance x'. 7. Procédé d'estimation de coefficients aérodynamiques d'un élément aérodynamique d'un aéronef dans un écoulement, ledit élément aérodynamique présentant, entre le bord d'attaque et le bord de fuite de celui-ci, un angle de flèche local dont la valeur varie entre une valeur de bord d'attaque et une valeur de bord de fuite, caractérisé en ce qu'il comprend les étapes de : - mise en oeuvre du procédé de simulation par 30 ordinateur du comportement aérodynamique dudit élémentaérodynamique selon l'une quelconque des revendications 1 à 6 ; puis - estimation des coefficients aérodynamiques dudit aéronef à partir de ladite solution numérique obtenue. 5 8. Procédé d'estimation selon la revendication 7, caractérisé en ce que ladite étape d'estimation des coefficients aérodynamiques comprend une étape intermédiaire d'estimation des coefficients locaux de 10 pression et de frottement pariétal à partir de ladite solution numérique obtenue et de la pression dynamique définie à partir de la vitesse amont incidente dudit écoulement. 15 9. Procédé d'estimation selon la revendication 8, caractérisé en ce que lesdits coefficients aérodynamiques locaux sont intégrés suivant ledit profil initial à l'épaisseur non modifiée, de manière à obtenir des coefficients aérodynamiques globaux 20 exprimés dans le repère lié à l'aéronef. 10. Procédé d'estimation selon la revendication 9, caractérisé en ce que lesdits coefficients aérodynamiques globaux ainsi obtenus sont projetés dans 25 un référentiel aérodynamique orienté suivant ladite vitesse amont, de manière à obtenir les coefficients de portance, de traînée et du moment de tangage. 30