These Jourieh

Développement d’un modèle représentatif d’une éolienne afin d’étudier l’implantation de plusieurs machines sur un parc éolien Munif Jourieh To cite this version: Munif Jourieh. Développement d’un modèle représentatif d’une éolienne afin d’étudier l’implantation de plusieurs machines sur un parc éolien. Sciences de l’ingénieur [physics]. Arts et Métiers ParisTech, 2007. Français. <NNT : 2007ENAM0032>. <pastel-00003390> HAL Id: pastel-00003390 https://pastel.archives-ouvertes.fr/pastel-00003390 Submitted on 8 Feb 2008 HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. N˚ : 2007 ENAM 0032 École doctorale n˚432 : Science des Métiers de l’Ingénieur THÈSE pour obtenir le grade de Docteur de l’École Nationale Supérieure d’Arts et Métiers Spécialité : MÉCANIQUE présentée et soutenue publiquement par Munif JOURIEH Développement d’un modèle représentatif d’une éolienne afin d’étudier l’implantation de plusieurs machines sur un parc éolien. à soutenir le 20 décembre 2007 à L’ENSAM, devant le Jury composé de : F. BAKIR Professeur, LEMFI, ENSAM, Paris Président G. DESCOMBES Professeur, CNAM, Paris Rapporteur I. AMJEH Professeur, Université de Damas Rapporteur F. MASSOUH Maı̂tre de Conférences HDR, LMF, ENSAM, Paris Examinateur P. KUSZLA Professeur agrégé, SINUMEF, ENSAM, Paris Examinateur G. NOTTON Maı̂tre de Conférences HDR, Université de Corse, Ajaccio Examinateur M. RAPIN Ingénieur de recherche, ONERA, Chatillon Invité Laboratoire de Simulation Numérique en Mécaniques des Fluides ENSAM CER de PARIS L’ENSAM est un Grand Établissement dépendant du Ministère de l’Éducation Nationale, composé de huit centres : AIX-EN-PROVENCE ANGERS BORDEAUX CHÂLONS-EN-CHAMPAGNE CLUNY LILLE METZ PARIS Table des matières 1 Introduction générale 1.1 Evolution de l’énergie éolienne . . . . . . . . . . . . . . 1.2 Sillage éolien et interaction entre machines . . . . . . . 1.3 Prédiction du comportement d’un rotor aérodynamique 1.4 Modèles hybrides . . . . . . . . . . . . . . . . . . . . . 1.5 Conclusion : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . : état de . . . . . . . . . . 2 Exploration en soufflerie du champ de vitesse autour d’une 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Description des moyens d’essais . . . . . . . . . . . . . . . . . 2.2.1 Soufflerie . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 Éolienne testée . . . . . . . . . . . . . . . . . . . . . . 2.3 Mesures PIV (Particle Image Velocimetry) . . . . . . . . . . . 2.3.1 Mise en place des mesures PIV . . . . . . . . . . . . . 2.3.2 Résultats . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.3 Champs de vitesse . . . . . . . . . . . . . . . . . . . . 2.3.4 Analyse des résultats . . . . . . . . . . . . . . . . . . . 2.4 Mesures par fil chaud . . . . . . . . . . . . . . . . . . . . . . 2.4.1 Principe de l’anémométrie à fil chaud . . . . . . . . . . 2.4.2 Méthodologie suivie dans l’exploration au fil chaud . . 2.4.3 Plan d’exploration . . . . . . . . . . . . . . . . . . . . 2.4.4 Représentation du sillage . . . . . . . . . . . . . . . . . 2.4.5 Étude de trajectoire des tourbillons marginaux . . . . . 2.4.6 Effet de l’éolienne sur l’écoulement à l’amont du rotor . 2.4.7 Champs de vitesse tangentielle . . . . . . . . . . . . . . 2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 . . . . . . l’art . . . . . . . . . . . . . . . . éolienne . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 8 12 14 17 18 . . . . . . . . . . . . . . . . . . 21 22 22 22 24 25 27 29 30 32 35 35 37 37 39 40 44 46 48 3 Développement d’un modèle hybride basé sur le concept du disque actif 3.1 Introduction : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Méthode de l’élément de pale . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Théorie de l’élément de pale . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Bilan de quantité de mouvement appliqué à un élément de pale . . 3.3 Application de la méthode de l’élément de pale . . . . . . . . . . . . . . . . 3.3.1 Algorithme de calcul . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Discrétisation de la pale et critère de convergence . . . . . . . . . . 3.4 Modèle de disque actif . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Equation de continuité . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2 Bilan de quantité de mouvement . . . . . . . . . . . . . . . . . . . 3.4.3 Coefficient de puissance . . . . . . . . . . . . . . . . . . . . . . . . 3.4.4 Limite de Betz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Description du modèle hybride . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Etude de la discrétisation du disque actif . . . . . . . . . . . . . . . . . . . 3.7 Interfaçage avec le logiciel Fluent . . . . . . . . . . . . . . . . . . . . . . . 3.8 Répartition volumique des forces appliquées sur le disque actif . . . . . . . 3.9 Algorithme de calcul . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Application du modèle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 Résultats et analyse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11.1 Eoliennes comparées . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11.2 Le domaine de calcul . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11.3 Comparaison des résultats numériques et expérimentaux . . . . . . 3.12 Comparaison entre le sillage mesuré d’une éolienne et le sillage calculé par le modèle hybride . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Développement d’un modèle hybride basé sur le concept de active 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Développement d’un modèle de cylindre actif . . . . . . . . . . . . 4.3 Répartition des efforts dans un cylindre actif . . . . . . . . . . . . 4.3.1 Choix du diamètre du cylindre actif . . . . . . . . . . . . . 4.4 Application du modèle de cylindre actif aux cas NREL s809 phase 4.4.1 Algorithme de calcul . . . . . . . . . . . . . . . . . . . . . 4.4.2 Surface de référence . . . . . . . . . . . . . . . . . . . . . . 3 51 52 52 53 55 57 58 58 62 62 63 63 64 64 66 68 68 70 73 73 73 75 77 78 80 la ligne . . . . . . . . . . . . . . . . . . . . II et VI . . . . . . . . . . 83 84 84 87 88 96 96 99 4.5 4.6 4.7 4.4.3 Domaine de calcul . . . . . . . . . . . . . . . . . . . . . . . . . . Résultats et analyse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Comparaison simulation/ expérience du sillage d’une éolienne Rutland 503 modifiée . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Interactions entre des éoliennes 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Simulation complète de l’interaction entre deux éoliennes à l’aide modèle hybride . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Effet de l’éloignement . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Effet de la vitesse à l’infini amont . . . . . . . . . . . . . . . . . . . 5.4.1 Effet de la vitesse en amont avec un modèle hybride CA . . 5.4.2 Effet de la vitesse en amont avec un modèle hybride DA . . 5.5 Interaction entre six machines en cascade . . . . . . . . . . . . . . . 5.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 . 101 . 103 . 105 108 . . . . 109 d’un . . . . 110 . . . . 112 . . . . 118 . . . . 118 . . . . 121 . . . . 122 . . . . 123 6 Conclusion générale 125 7 Annexe 130 4 Liste des symboles V1 : Vitesse du vent à l’infini amont. [m/s] Vd : Vitesse du vent dans le plan du rotor . [m/s] V2 : Vitesse du vent à l’infini aval. [m/s] Vl : vitesse locale. [m/s] W : Vitesse relative du vent. [m/s] D : Diamètre du rotor. [m] d : Diamètre du cylindre actif (CA). [m] R : Rayon de l’extrémité de la pale. [m] r : Rayon local. [m] B : Nombre des pales. [−] Cx : Coefficient de traı̂née. [−] Cz : Coefficient de portance. [−] Cn : Coefficient de la force normale au plan de rotation. [−] Ct : Coefficient de la force tangentielle au plan de rotation. [−] Cp : Coefficient de puissance. [−] β : Angle de calage. [rad] α : Angle d’incidence. [rad] φ : Angle d’écoulement relatif. [rad] θ : Angle azimutal. [rad] ρ : Masse volumique. [kg/m3 ] ν : Viscosité cinématique. [m2 /s] Ω : Vitesse angulaire du rotor. [rad/s] N : Vitesse de rotation du rotor. [tr/min] c : corde de pale.[m] a : Facteur d’induction axial. [−] a′ : Facteur d’induction tangentiel. [−] P + : Pression en amont du disque actif. [P a] P − : Pression en aval du disque actif. [P a] P1 : Pression atmosphérique. [P a] A1 : Surface du tube de courant à l’infini amont. [m2 ] Ad : Surface du disque actif. [m2 ] A2 : Surface du tube de courant à l’infini aval. [m2 ] σ : Solidité de la pale. [−] λ : Rapidité spécifique. [−] 5 P : Puissance extraite par une éolienne. [W ] 6 Chapitre 1 Introduction générale 7 1.1 Evolution de l’énergie éolienne L’énergie éolienne est probablement une des plus anciennes sources d’énergie. Cette énergie propre et renouvelable existe depuis toujours, mais jusqu’à présent son exploitation reste difficile. L’utilisation de l’énergie éolienne a commencé en 1700 avant Jésus-Christ environ, Hammourabi, roi de Babylone, actuellement l’Irak, a pensé bénéficier de cette énergie propre pour l’irrigation. Il a utilisé la puissance du vent pour le pompage de l’eau avec des éoliennes à axe vertical [1, 2]. Au neuvième siècle avant Jésus-Christ, la Perse, royaume qui correspond maintenant à l’est de l’Iran et à l’Afghanistan, a utilisé l’énergie du vent dans le but de moudre le grain et de pomper l’eau. La première conception documentée connue était celle d’un moulin à vent persan, le panémone, utilisant des voiles verticales faites de roseaux ou de bois qui étaient attachées à un axe central vertical. La figure 1.1 montre le principe de fonctionnement d’un tel moulin [3, 4]. Vent Voile Fig. 1.1 – Principe du premier moulin à vent persan à axe vertical. En Inde, l’exploitation de cette énergie renouvelable a pris son essor au quatrième siècle avant Jésus-Christ [5]. 8 Ensuite, trois siècles avant Jésus-Christ, les Egyptiens ont commencé à bénéficier de cette énergie propre. L’inventeur égyptien Héron d’Alexandrie utilisa l’énergie éolienne grâce à un moulin à vent à axe horizontal. Ce dispositif permettait, contrairement à ceux utilisés pour les tâches agricoles, d’utiliser l’énergie éolienne pour alimenter en air comprimé un orgue [1]. Au Moyen-Age, en Europe, l’exploitation de cette énergie a commencé par l’apparition des moulins de vent, principalement en France, en Italie, en Espagne et au Portugal. Les premiers travaux écrits ont ainsi vu le jour en Normandie en 1180. Certains auteurs pensent que la technologie du moulin à vent à axe vertical développée en Perse s’étendit en Europe suite aux Croisades au Moyen-Orient en évoluant sous la forme de moulins à vent à axe horizontal. On les rencontre, un peu plus tard, en Grande-Bretagne, en Hollande et en Allemagne sous le type de machine à axe horizontal comportant quatre ailes placées en croix. Elles servaient principalement à moudre du grain [6]. Les Hollandais participèrent activement au développement des moulins en Europe, grâce aux nombreuses améliorations dans la conception et à l’invention de différents types de moulins. En effet, pendant le treizième siècle, la Hollande utilisa les moulins à vent pour pomper l’eau et ainsi assécher les polders. Accouplées à des roues à godet ou à des vis d’Archimède, ces machines pouvaient élever l’eau jusqu’à cinq mètres. L’utilisation des moulins à vent a connu un grand succès jusqu’au dix neuvième siècle. Ils ont fourni à l’homme l’énergie mécanique nécessaire pour certains travaux agricoles. Mais avec la révolution industrielle et l’apparition de la machine à vapeur, du moteur à combustion et plus tard le développement de l’électricité, le développement des éoliennes et l’exploitation des moulins à vent sont délaissés. Leurs utilisations ont décliné jusqu’après la deuxième guerre mondiale. Face à la présence de nouveaux moyens de production d’énergie, les aérogénérateurs n’arrivent pas à s’imposer [2]. Pendant le vingtième siècle, plusieurs projets d’éoliennes voient le jour. L’éolienne lente multipale est développée en Amérique par la Rural Electrification Administration [7]. L’éolienne rapide, inventée en France par l’Académicien français Darrieus, entraı̂nait des générateurs électriques [1]. En 1950, Johannes Juul [6] développa un modèle éolien avec trois pales, utilisant des dispositifs de réglages aérodynamiques de la puissance dans le cas de décrochage et du contrôle de dérapage. Cependant dans les années 1970 et après la crise du pétrole causée par la guerre au Moyen-Orient, l’administration du président Carter aux USA et d’autres dans plusieurs pays comme la Suède, le Canada et la Grande-Bretagne, ont pris la décision de développer la recherche dans le domaine de l’énergie renouvelable [8, 9]. En effet, avec la diminution du stock mondial d’hydrocarbures, la demande énergétique sans cesse croissante, et la crainte 9 d’une pollution de plus en plus envahissante, l’énergie propre et renouvelable attire les puissances mondiales. L’énergie éolienne revient au premier plan de l’actualité car son exploitation peut s’avérer très rentable dans les nombreuses régions ventées du globe. Ainsi, au début des années 1980, dans des pays tels que le Danemark ou l’Allemagne, où il n’existe pas des ressources d’énergie aussi importantes que l’énergie éolienne, le développement d’une industrie éolienne performante a été privilégié [10, 11, 7]. Les éoliennes ont ainsi continué à évoluer au cours des 20 dernières années, et le coût global de l’énergie nécessaire à la production d’électricité à partir du vent est maintenant concurrentiel avec les sources d’énergie traditionnelles comme les combustibles fossiles [12, 13]. Cette réduction du coût de l’électricité est le résultat de progrès importants de la technologie utilisée par cette industrie (amélioration des conceptions aérodynamiques, amélioration des matériaux utilisés, . . .) [14, 15]. Actuellement, l’énergie éolienne est bien implantée parmi les autres sources d’énergie avec une croissance très forte [16, 17]. La figure 1.2 présente la puissance issue de l’énergie éolienne des dix premiers marchés éoliens du monde en 2005 [18]. A partir de celle-ci, on observe que cette énergie est bien développée aux USA, Allemagne, Espagne et en Inde par rapport aux autres pays. Fig. 1.2 – Production d’énergie éolienne dans les dix premiers marchés éoliens en 2005. 10 Fig. 1.3 – Puissance et croissance du secteur de l’énergie en France (MW, MW/an). La croissance de l’exploitation de l’énergie éolienne dans le monde est influencée par un facteur très important, c’est la pollution, car les autres types d’énergie (Pétrole, Nucléaire) sont forcément polluants, ce qui amène à prendre des initiatives dans le monde pour ajouter des conditions d’utilisation afin de diminuer cette pollution (protocole de Kyoto [19, 20]). Ce qui donne des avantages supplémentaires à l’utilisation de l’énergie éolienne. La figure 1.3 représente l’évolution de l’énergie éolienne en France de 1996 à 2005. On remarque une forte croissance de l’exploitation de l’énergie éolienne au cours des cinq dernières années. Étant donnée cette forte croissance et le fait que les zones où l’énergie éolienne est techniquement exploitable sont limités, il est important de bien gérer les parcs éoliens. L’installation d’un parc éolien doit donc être optimisée de manière à tirer parti au mieux de la ressource en vent. Cette recherche de l’installation optimale fait l’objet de cette thèse où nous allons tenter de mettre en place des modèles de rotors d’éoliennes qui permettront de représenter un parc de machines et d’en évaluer les performances. 11 1.2 Sillage éolien et interaction entre machines Dans les parcs éoliens, il est nécessaire d’espacer les éoliennes afin d’éviter que le sillage et le déficit de vitesse existant derrière chaque machine n’affectent trop la production énergétique et l’intégrité mécanique des éoliennes situées plus en aval [21, 22]. En règle générale, les éoliennes dans les parcs sont toutes identiques et la distance entre les machines est définie d’après des règles simples données par les constructeurs, à savoir : cinq à neuf fois le diamètre dans la direction des vents dominants et de trois à cinq fois le diamètre dans la direction perpendiculaire [23]. Ces règles ne constituent évidemment pas une garantie d’optimisation de l’implantation des machines et ne permettent pas d’évaluer la prise en compte de la topographie du terrain. Actuellement, la simulation numérique en mécanique des fluides a fait d’énormes progrès et permet de donner des informations intéressantes sur l’écoulement d’air autour des éoliennes. Les codes industriels modernes permettent d’obtenir des résultats quantitatifs satisfaisants [24]. La simulation pourrait donc être un outil appréciable pour analyser le sillage éolien [25] et optimiser l’implantation des éoliennes dans un parc en fonction des caractéristiques des machines et de la topographie du site éolien [26, 27, 28]. Cependant, la simulation de l’écoulement à travers le rotor éolien et dans son sillage est toujours une opération lourde nécessitant une préparation et un temps de calcul important. Les moyens de calcul modernes permettent de simuler finement le comportement aérodynamique d’une pale d’éolienne et éventuellement d’un rotor complet [29]. Il est toutefois hors de question, dans l’état actuel du matériel informatique disponible, de simuler l’ensemble des machines composant un parc éolien. D’un point de vue expérimental, la non uniformité du champ de vitesse sur un site éolien, de même que l’encombrement des machines rend pratiquement impossible l’étude en soufflerie de l’ensemble d’un parc éolien. Ainsi, la voie expérimentale se limite à des études partielles des parcs éoliens avec un nombre très réduit de machines. Il s’agit donc dans le travail réalisé dans cette thèse de proposer un modèle capable de dépasser les limitations imposées tant par le dispositif expérimental que par les outils de simulation. Cette thèse partira de l’idée qui consiste à rechercher un modèle équivalent du rotor éolien permettant de bien représenter le sillage et son développement tout en évitant de calculer en détail l’écoulement autour des différentes pales. Le temps de calcul ainsi épargné, grâce à l’usage d’un tel modèle, doit permettre de simuler et d’optimiser l’implantation des machines d’un parc éolien dans un environnement donné. L’approche habituellement utilisée pour modéliser les rotors éoliens et les rotors d’hélicoptères aussi, est basée sur un disque actif [30, 31, 32]. Le rotor complet est ainsi remplacé par un 12 disque infiniment mince à nombre infini de pales et qui génère un écoulement équivalent à l’écoulement réel moyenné. Concrètement, c’est une condition aux limites particulière (discontinuité de pression) appliquée sur toute la surface balayée par les pales, condition qui produit la même variation de quantité de mouvement et d’énergie que le rotor en fonctionnement. L’objectif recherché consiste à trouver un meilleur modèle que le disque actif et qui rendra possible la prise en compte du nombre fini des pales du rotor. Ainsi, il est question de développer une modélisation capable de reproduire chaque pale individuellement. 13 1.3 Prédiction du comportement d’un rotor aérodynamique : état de l’art Nous allons ici faire un bref tour d’horizon des méthodes utilisées dans la simulation des écoulements autour des rotors éoliens. La plupart des méthodes appliquées pour prévoir le comportement d’une éolienne dans le cas axisymétrique, utilisent la théorie du disque actif, où le rotor est remplacé par des forces distribuées sur un disque d’épaisseur nulle et perméable [33, 4]. Cette théorie a été introduite par Froude [34] et Rankine [35], elle est développée ensuite par Lanchester en 1915 [36] puis par Betz en 1920 [5]. Ils ont démontré 16 que la puissance maximale extraite par une éolienne est égale à 27 de l’énergie cinétique de l’écoulement incident, soit environ 59.3%. En 1935, la théorie de l’élément de pale est proposée par Glauert [37], cette méthode est basée sur la division de l’écoulement en volumes de contrôles annulaires, auxquels on applique le bilan de quantité de mouvement et d’énergie. Ces anneaux s’étendent de l’infini amont jusqu’à l’infini aval par rapport au rotor, figure 1.4. Y Disque actif V1 R r Z X dr Fig. 1.4 – Tube de courant au rayon local r à travers une éolienne tripales remplacée par un disque actif. Du fait de sa simplicité, cette méthode est largement utilisée dans le domaine industriel [38]. L’application de la théorie de l’élément de pale sur le disque actif donne la base la plus utilisable pour construire des modèles hybrides [39, 40, 41]. 14 La limitation principale de la méthode du disque actif est qu’elle n’est valide que pour les cas des éoliennes axisymétriques, et que cette méthode distribue les forces de façon égale dans la direction azimutale du disque actif et ne représente pas individuellement chaque pale. Ainsi, l’influence des pales est prise comme une quantité intégrée dans la direction azimutale. Pour surmonter ce problème et prendre en compte la présence des pales, Sørensen et Shen [42], ont introduit la technique de la ligne active, où les forces créées par une pale sont distribuées le long d’une ligne qui représente cette pale, figure 1.5. Dans le concept de la ligne active, chaque pale est représentée par une ligne discrétisée par n points. Les forces sont distribuées autour de chaque point dans un plan perpendiculaire à la ligne active selon une distribution spécifique. Par exemple, une distribution Gaussienne dans le travail de Mikkelesen [43], ou une distribution uniforme dans le modèle de cylindre actif présenté dans cette thèse au chapitre quatre. y x Ligne active Fig. 1.5 – Concept de ligne active La première formulation du concept de ligne active a été proposée par Sørensen et Shen [42], leur analyse montre l’accord avec la courbe de puissance mesurée pour l’éolienne 15 Nordtank tri pale 500kW. Plus tard, Mikkelsen [43] a proposé un modèle hybride en combinant cette méthode avec le code EllipSys-3D. Les résultats obtenus sont comparés avec les résultats expérimentaux de la turbine de Tjaereborg 1 pour plusieurs cas de dérapage, et donnent des valeurs assez proches et comparables [39]. Ivanell [44] a utilisé la même méthode que Mikkelsen et Sørensen, et a montré que le contrôle de la distribution Gaussienne des forces autour de la ligne active est influencé par quelques paramètres. Il a mis en évidence l’influence du nombre de Reynolds et de la résolution du maillage. Ivanell a aussi observé l’amélioration des résultats en utilisant la correction de Prandtl afin de mieux prendre en compte les phénomènes d’extrémité de pale. D’autre part, Massouh et Dobrev [45] ont proposé un nouveau modèle hybride, dans lequel chaque pale est remplacée par une surface de discontinuité de pression. De cette manière, on surmonte les difficultés de la distribution des forces autour de la ligne active. Cette distribution de pression varie le long de la corde. La représentation de la pale et des conditions initiales pour le développement du sillage est ainsi améliorée par rapport aux méthodes de disque actif ou de ligne active. La validation de ce modèle hybride est faite avec le cas de l’éolienne NREL2 s8093 phase IV et aboutit à des résultats satisfaisants. On observe que les méthodes basées sur le concept du disque actif ou sur le concept de la ligne active et la théorie de l’élément de pale BEM (Blade Element Momentum), ont bien validé leur capacité à représenter différents types de rotors(pales droites ou vrillées) dans différentes conditions de fonctionnement (dérapage par exemple) [46, 47, 48]. Toutes les méthodes trouvées dans la littérature sont des représentations plus ou moins approchées des rotors réels. De ce fait on observe une dispersion importante des résultats fournis par chacun des auteurs. Ceci est particulièrement visible avec la comparaison simulation-expérience établie par le NREL en Décembre 2000 [49]. Dans cette étude, 19 méthodes de simulation numérique d’une éolienne NREL phase VI sont comparées. Elles utilisent différents modèles, de la méthode la plus simple, la méthode de l’élément de pale (BEM) jusqu’aux simulations tridimensionnelles turbulentes complètes. La figure 1.6, représente les modèles participants à la comparaison pour une vitesse de rotation faible et sans dérapage, la ligne noire avec les points représente les mesures expérimentales pour 1 La turbine à axe horizontal de Tjaereborg, d’un diamètre de 61 m, a trois pales. Le profil des pales est un NACA 4412-43. La longueur de la corde vaut 0.9 m à l’extrémité de pale et augmente jusqu’à 3.3 m au rayon extérieur de 6 m. L’axe de rotation est situé à 60 m du sol. Le vrillage de pale est de 1° tous les 3 m. 2 NREL : National Renewable Energy Laboratory. 3 Pour plus des information autour le profil s809, voir l’annexe A. 16 l’éolienne NREL phase VI dans la soufflerie Ames de la NASA4 [50]. On voit nettement sur cette figure qu’il y a plusieurs méthodes qui ont donné des résultats assez loin des données expérimentales. Même d’un simple point de vue qualitatif, certaines approches semblent mises en echec. Ces résultats mettent en évidence la difficulté de produire des simulation prédictives. Fig. 1.6 – Présentation des méthodes de prédiction de couple de l’éolienne NREL phase VI en fonction de la vitesse du vent. 1.4 Modèles hybrides Notre contribution à la simulation des écoulements autour des éoliennes consistera à développer, dans la lignée des méthodes existantes un nouveau modèle hybride. C’est à dire que nous tenterons de combiner des modèles de rotor (éléments de pale, disque actif, cylindre actif) avec une simulation numérique basée sur la résolution des équations de Navier-Stokes tridimensionnelles turbulentes. Dans ce travail nous présentons deux types de modèles hybrides, le premier est basé sur le concept du Disque Actif (DA), et le deuxième est basé sur le concept de la ligne active et sera un modèle de Cylindre Actif (CA). 4 NASA : National Aeronautics and Spase Administration. 17 – Avec le modèle hybride DA, nous représentons le rotor par un disque de diamètre égal à celui du rotor et d’épaisseur liée à la corde de la pale5 . La distribution des forces exercées par les pales est faite de manière uniforme sur le disque. La présence des pales est donc moyennée. Ce type de modèle hybride n’est utilisable que dans les cas axisymétriques et ne peut pas calculer un rotor en dérapage. – Pour surmonter cette limitation, nous utilisons ensuite un modèle hybride dit modèle de cylindre actif (CA) qui prend en compte le nombre de pales. Ce modèle est basé sur le concept de la ligne active, chaque pale est remplacée par un cylindre de longueur égale à la longueur de la pale et de diamètre lié à la corde de la pale. La validation de ces modèles hybrides est faite de deux manières différentes, par rapport au sillage éolien et par rapport à l’évolution de la puissance en fonction de la vitesse du vent. Pour le sillage éolien, nous avons réalisé des travaux expérimentaux pour une éolienne tri pales (Rutland 503 modifiée) dans la soufflerie du laboratoire de mécanique des fluides (LMF) à l’École Nationale Supérieure d’Arts et Métiers (ENSAM) à Paris. Ce travail présente l’exploration du sillage éolien dans la veine d’essai de la soufflerie, celui-ci est exploré par deux méthodes différentes. La première méthode utilisée est l’anémométrie à fil chaud à deux composantes, cela donne l’exploration du champ des vitesses axiale et tangentielle du sillage. La deuxième méthode est la méthode PIV (Particle Image Velocity). Avec celle-ci, on explore le champ des vitesses axiales et radiales. Pour avoir une présentation assez détaillée sur le sillage derrière l’éolienne, on a fait des explorations à plusieurs angles azimutaux. Pour l’évolution de puissance en fonction de la vitesse du vent, la validation de ces modèles est faite avec une éolienne NREL phase II [51] et NREL phase VI [52]. Enfin, nous avons étudié l’interaction entre des machines en utilisant ces modèles hybrides. 1.5 Conclusion : Ce chapitre donne une idée générale sur de l’énergie éolienne et son intérêt. Cette énergie renouvelable a suivi son chemin depuis plusieurs années avec une croissance annuelle d’utilisation très importante dans le monde. Au début du vingtième siècle, plusieurs hypothèses sont apparues dans le domaine de l’aérodynamique concernant les éoliennes. Les théories du disque actif et la ligne active sont les théories les plus utilisées pour étudier 5 Par abus de langage, nous retiendrons ici l’expression ”disque actif” bien que ce disque possède une épaisseur non nulle et soit en fait un cylindre de hauteur très petite devant son rayon. 18 la performance d’une éolienne avec les équations de Navier-Stokes. La progression dans le domaine de la simulation numérique permet de représenter une géométrie reélle d’une éolienne par un modèle hybride beaucoup plus simple et facile à réaliser. 19 20 Chapitre 2 Exploration en soufflerie du champ de vitesse autour d’une éolienne 21 2.1 Introduction Cette partie présente des explorations de l’écoulement en amont et en aval d’une éolienne tri pales à axe horizontal. Ces travaux permettent d’avoir les données expérimentales nécessaires pour la validation de nos modèles hybrides, validation qui sera présentée dans les chapitres suivants. Dans ce travail, nous avons utilisé deux méthodes différentes, la première est la méthode de vélocimétrie par image de particules PIV (Particle Image Velocimetry). Par cette méthode, nous avons exploré le sillage proche derrière la machine testée. Les mesures ont été synchronisées avec la position d’une pale de référence. Pour avoir une idée claire de la structure du sillage en trois dimensions, d’autres mesures ont été réalisées avec un retard contrôlé afin d’explorer plusieurs plans azimutaux derrières la pale, l’angle azimutal entre deux plans successif est égal à 30°. Avec cette méthode, nous n’avons pas pu explorer les deux champs de vitesse axiale et tangentielle simultanément. Pour atteindre cet objectif, nous avons utilisé une autre méthode, la méthode de l’anémométrie à fil chaud. Avec l’anémométrie à fil chaud, nous avons fait une exploration pour les deux champs de vitesse (axiale et tangentielle). Cette exploration est réalisée en amont du plan de rotation jusqu’à 0.4D1 , et en aval jusqu’à 1.5D environ du celui-ci. Ces expériences sont réalisées dans le veine d’essai de la soufflerie de l’ENSAM-Paris. A cause des dimensions limitées de notre veine d’essai, il est possible d’explorer le sillage proche uniquement. Mais, l’intérêt ici est d’avoir des données expérimentales concernant le sillage proche, et que ces données permettent de valider des modèles de rotors simplifiés, couplés avec la simulation numérique. De plus, le sillage proche donne les conditions initiales pour le développement du sillage lointain [53, 54]. 2.2 Description des moyens d’essais 2.2.1 Soufflerie La soufflerie du laboratoire de mécanique des fluides de l’ENSAM-Paris est de type Prandtl avec une veine semi guidée sans parois latérales et ouverte à la pression atmosphérique. Elle comporte les parties suivantes : veine d’essai, diffuseur, deux séries de coudes avant la veine de retour, ventilateur, et enfin deux séries de coudes avant la chambre de tranquillisation puis le convergent d’alimentation de la veine, voir la figure 2.1. 1 D : Diamètre du rotor testé. 22 La veine d’essai a une section rectangulaire 1.35 m × 1.65 m avec une longueur de 2 m. Cette partie est ouverte à la pression atmosphérique et la vitesse maximale de l’écoulement est 40 m/s. La veine d’essai est dotée d’un plafond transparent en plexiglas qui facilite la visualisation des écoulements. L’équipement de la veine d’essai comporte un explorateur commandé par ordinateur et capable de se déplacer dans les trois directions de l’espace. La veine est suivie par un diffuseur avec un angle d’ouverture voisin de 7 degrés pour éviter les décollements de l’écoulement sur les parois. Chambre de tranquillisation. Veine d’essai 1.35m x 1.65m x 2m. 7 m. 8.4 m. Explorateur. 22.5m. Ventilateur électrique 75kW. 3 pales. Explorateur. Veine de retour 3m x 3m x 6m. Fig. 2.1 – Soufflerie de l’ENSAM-Paris. La veine de retour possède une section carrée de 3 m × 3 m sur une longueur de 6 m et la vitesse maximale d’écoulement est d’environ 10 m/s. Cette veine est intéressante par ses dimensions, mais elle n’offre pas de bonnes caractéristiques d’uniformité du champ de vitesses et pour ces raisons, l’exploration du champ de vitesse dans la veine d’essai est préférable. Le ventilateur comporte une hélice tripale de trois mètres de diamètre pouvant fonctionner jusqu’à 750 tour/min. Il est entraı̂né par un moteur électrique de 75 kW . La variation de vitesse de rotation du moteur permet d’assurer la régulation de la vitesse d’écoulement. Après le ventilateur, l’air est orienté vers la chambre de tranquillisation qui comporte une grille permettant d’aligner les filets fluides et de diminuer la turbulence. La chambre de tranquillisation alimente la veine d’essai par l’intermédiaire d’un convergent qui permet une répartition uniforme des vitesses dans sa section contractée. Le convergent 23 entre la chambre de tranquillisation et la veine d’essai possède un rapport de contraction de 1 :12.5, ce qui permet de rendre l’écoulement globalement uniforme. Ainsi, le vent atteint dans la veine d’essai une vitesse maximale de 40 m/s (144 km/h) avec une très bonne qualité d’écoulement. Enfin, il y a quatre coudes à 90° munis d’aubages qui imposent au fluide un écoulement en plans parallèles. Les coudes ont pour effet de guider au mieux l’écoulement pour éviter la formation de structures tourbillonnaires lors des changements de direction. 2.2.2 Éolienne testée La maquette utilisée dans ce travail est une éolienne commerciale (Rutland 503 de la Société Marlec) qui a été modifiée pour les essais, figure 2.2. En effet, l’éolienne d’origine avait 6 pales, mais trois pales ont été ôtées afin d’obtenir une éolienne tripales semblable aux éoliennes utilisées dans les fermes éoliennes. Le diamètre de cette éolienne est de 0.5 m avec un moyeu de 0.135 m de diamètre. Les pales ne sont pas vrillées et sont montées avec un angle de calage constant de 10°. La corde au pied de la pale est 0.065 m et 0.045 m à l’extrémité. L’axe de rotation de cette éolienne est placé à une hauteur de 0.7 m au milieu de la veine d’essai, à l’aide d’un mat de diamètre 0.037 m pour éviter l’influence du plafond et du plancher sur le comportement du sillage et pour permettre aux lasers fixés sur le plafond transparent de la veine d’essais d’éclairer le plan d’exploration avec une intensité suffisante. L’exploration du sillage de cette maquette est réalisée avec une vitesse axiale de 9.4 m/s et une vitesse de rotation de 110 rad/sec. Le contrôle de la vitesse de rotation a été assuré par un rhéostat connecté à la sortie de la génératrice électrique équipant l’éolienne. 24 Fig. 2.2 – L’éolienne testée dans la veine d’essai. 2.3 Mesures PIV (Particle Image Velocimetry) Depuis une vingtaine d’années, la technique PIV est utilisée pour mesurer la vitesse d’un fluide et déterminer le champs de vorticité instantanée. Cette technique permet de mesurer le champ instantané de vitesse au lieu de mesurer la vitesse en un point en fonction du temps. Le principe de cette méthode est basé sur l’ensemenssement de l’écoulement et la mesure, en fonction du temps, de la vitesse des particules et sur la détermination de leurs trajectoires dans un plan illuminé par une nappe laser [55]. L’illumination intensive par le laser d’un plan de l’écoulement, permet d’avoir une image instantanée, et de détecter les positions instantanées des particules d’ensemencement. On enregistre les images des particules à deux instants successifs. La synchronisation entre les tirs laser et la prise d’images permet d’obtenir un champ de vitesse instantanée. Le principe de mesure est expliqué sur la figure 2.4 qui représente un plan rectangulaire détecté par une caméra. 25 Fig. 2.3 – Éolienne testée avec une nappe Laser. Pour simplifier, on suppose la présence de 6 particules qui passeront des positions indiquées A,B,C,D,E et F dans la première image aux positions A1,B1,C1,D1,E1 et F1 dans la deuxième image pendant un temps égal à △t. B1 A1 A B D1 C1 D C E1 E F F1 Fig. 2.4 – Déplacement des particules. A partir de cette figure 2.4, les deux composantes de la vitesse instantanée selon X et Y sont calculables. VX = △X △ t. VY = △Y  △ t. △t : Temps entre deux images successives (Le temps entre deux tirs laser). 26 △X : Distance entre la première et la deuxième position pour une particule selon l’axe X. △Y : Distance entre la première et la deuxième position pour une particule selon l’axe Y. V : Vitesse instantanée pour une particule dans la tranche éclairée. En réalité, l’image PIV contient un grand nombre de particules. Un calcul d’intercorrélation est utilisé pour identifier la position des particules dans les deux images successives. La répétition des mesures permet d’obtenir le champ de vitesses moyennées [56, 57]. 2.3.1 Mise en place des mesures PIV Dans la soufflerie de l’ENSAM-Paris, la partie en plexiglas transparente dans le plafond de la veine d’essai a permis d’installer les lasers Nd-Yag au-dessus de la veine. Comme la veine est semi guidée et sans parois latérales, la caméra a été installée à coté de la veine pour assurer la prise d’images dans le sens normal au plan d’éclairage laser. Ainsi, le matériel PIV est entièrement à l’extérieur de la veine d’essais et ne perturbe pas l’écoulement, voir la figure 2.5. L’ensemencement est réalisé grâce à l’utilisation de gouttelettes d’huile d’olive. Pour avoir un écoulement avec une répartition homogène des particules d’ensemencement dans la tranche éclairée par le laser, qui aide à bien représenter les vecteurs du champ de vitesse, un générateur est utilisé pour pulvériser des gouttelettes d’huile dans le diffuseur de la soufflerie. Le temps nécessaire pour parcourir la veine de retour et la chambre de tranquilisation a permis une excellente répartition des gouttelettes dans l’écoulement en arrivant à la veine d’essai. Pour synchroniser les mesures PIV avec la position angulaire du rotor, un codeur optique situé sur le mat de l’éolienne a été utilisé pour viser une cible tournante collée sur le moyeu et émettre un signal de référence temporelle qui indique la position angulaire du rotor. Ce signal permet de contrôler le déclenchement des tirs laser et la prise d’images. Le temps entre deux tirs laser a été réglé à : △ T = 150 µs. Pour explorer le volume occupé par le sillage et les trajectoires des tourbillons marginaux, nous avons mesuré le champ de vitesse dans plusieurs plans azimutaux avec les angles 0, 30, 60 et 90º ; le zéro correspondant à la position verticale de pale, voir la figure 2.5. A cause de la limite de la capacité de notre matériel, il est impossible d’explorer le champ de vitesse derrière la machine par une seule image. Pour surmonter ce problème, 27 Fig. 2.5 – Banc d’éssai. nous avons divisé le champ de vitesses en six fenêtres (3 horizontales × 2 verticales), voir la figure 2.6. L’échelle des fenêtres et leurs positions par rapport au rotor ont été définies à l’aide de mires placées dans le plan de l’exploration après chaque série d’essais. Ainsi, quatre séries de 6 fois 95 paires d’images ont été réalisées afin d’obtenir les différents plans d’exploration azimutaux définis avec les angles 0, 30, 60 et 90º ; figure 2.5. En raison de la vitesse de rotation élevée de l’éolienne (autour de 1050 tour/min), la caméra prenait une paire d’images tous les deux tours. Pour chaque fenêtre, la prise d’images était répétée 95 fois afin de reproduire une séquence temporelle d’environ 12 secondes et d’améliorer le calcul de la vitesse moyenne. 28 h2 h1 m2 m1 h3 m3 Fig. 2.6 – Six fenêtres d’exploration h1, h2, h3 (haut), et m1, m2, m3 (bas) et champs des vitesses moyennées. 2.3.2 Résultats La figure 2.7 présente une image brute prise dans la fenêtre h1 (1600 × 1200 pixels) directement derrière le rotor en position verticale de la pale. Cette image montre la pale (à gauche) et les noyaux des tubes tourbillonnaires issus des extrémités des autres pales. Pale R=0.25 m Noyaux Fig. 2.7 – Image brute, fenêtre h1. 29 En raison des vitesses induites par les tubes tourbillonnaires et qui écartent les particules d’ensemencement du centre du noyau, ce dernier apparaı̂t sur l’image comme un cercle noir par manque de particules éclairées. Le traitement des images brutes permet de constituer les champs des vitesses instantanées et moyennées. Après avoir obtenu un nombre suffisant d’images par la technique PIV, le traitement de celles-ci aura lieu par un algorithme itératif multi-résolutions, Susset [58]. Nous présentons les résultats dans les parties suivantes. 2.3.3 Champs de vitesse Dans ce paragraphe, nous allons représenter les champs des vitesses instantanées et moyennées. La figure 2.8 montre le champ de vitesse instantanée correspondant à l’image brute dans la fenêtre h1, déjà présentée sur la figure 2.7. Ce champ montre clairement l’intersection des tubes tourbillonnaires émanant des extrémités des pales et leurs interactions avec le plan d’exploration. Par ailleurs, on peut voir facilement l’augmentation du diamètre du tube de courant comme résultat du ralentissement de l’écoulement créé par l’éolienne, comme le prévoit la théorie de Froude-Rankine [59]. Fig. 2.8 – Champ de la vitesse instantanée, fenêtre h1. 30 Pour obtenir le champ de vitesse moyennée pour chaque fenêtre, on a moyenné les champs instantanés résultant des séries d’images prises (95 images). La figure 2.9 montre le champ de la vitesse moyennée pour la fenêtre h1. Fig. 2.9 – Champ de la vitesse moyennée, fenêtre h1. Il est à signaler que l’examen des champs successifs dans le temps montre une fluctuation de la position des noyaux des tourbillons marginaux au delà d’une distance de l’ordre 0.75 D par rapport au plan de rotation. Cette fluctuation est due à la variation temporelle de la puissance du rotor. En effet, lors d’une augmentation de la puissance absorbée, la vitesse moyenne dans le sillage diminue conduisant ainsi à une diminution du pas des tourbillons marginaux [56]. En conséquence, le respect de l’équation de continuité fait que le rayon du tube de courant augmente et déplace les tourbillons vers l’extérieur. On peut noter que la fluctuation radiale des noyaux dépend faiblement de leurs positions axiales. Par contre, leur battement axial augmente rapidement, parce ce que ce battement est amplifié en fonction du nombre de pas parcourus. En conséquence du battement du noyau, une réduction de l’intensité du tourbillon apparaı̂t lors du calcul de la moyenne des champs instantanés. Ainsi, le champ de vitesse stationnaire (champ moyenné) n’est pas tout à fait représentatif du développement du sillage. A l’aide des mesures des champs de la vitesse axiale moyennée dans la fenêtre h1 et m1 dans les quatre positions azimutales, 0, 30, 60 et 90°, une représentation tridimensionnelle 31 est obtenue pour le sillage proche du rotor, voir la figure 2.10. Nous observons à l’extrémité de chaque pale un tube tourbillonnaire qui sort en trajectoire hélicoı̈dale. Ainsi, cette figure 2.10, donne une idée de la forme de sillage et le volume occupé par celui-ci. 0° 30° 60° 90° Fig. 2.10 – Représentation 3D du sillage proche. 2.3.4 Analyse des résultats L’analyse des vitesses dans le sillage du rotor, figure 2.11, montre que le moyeu est une source de perturbations et d’instationnarités aérodynamiques fortes. Ceci vient de la forme encombrante du moyeu qui comporte la génératrice. Par ailleurs, le faible angle de calage au pied de la pale conduit à un blocage (avec le moyeu) important dans le chemin de l’écoulement. Par conséquence, un ralentissement important aura lieu derrière le moyeu, qui donne une raison de l’existence d’un courant de retour. L’incertitude sur les mesures de vitesse au voisinage de la pale vient des difficultés spécifiques de la technique PIV pour explorer l’écoulement près des parois, par exemple, 32 Courant de retour Fig. 2.11 – Lignes de courant dans les plans h1 et m1, angle azimutal 0°. V=0 Fig. 2.12 – Champ de la vitesse axiale induite par les tubes tourbillonnaires. la réflexion du laser sur la surface de la pale, la coupe du trajet suivi par un particule au voisinage de la pale et la valeur élevée de la vitesse tangentielle à quelques millimètres de la pale, voir figure 2.9. Pour étudier l’influence des tourbillons marginaux, nous avons présenté le champ de la vitesse induite par les tubes tourbillonnaires (vitesse axiale locale moins la vitesse en amont), le résultat est représenté sur la figure 2.12. Ainsi, la zone tracée par une ligne en pointillés au niveau des noyaux, avec une vitesse nulle, divise l’écoulement en deux 33 parties : une partie interne où la vitesse axiale est ralentie par les tourbillons marginaux et une partie externe où la vitesse est accélérée par les tourbillons. Cette image peut rappeler à un certain point un modèle tourbillonnaire simple (modèle de Prandtl) où le rotor est représenté par un disque de tourbillons annulaires et le sillage par un cylindre avec des tourbillons annulaires espacés d’une distance égale au pas de l’hélicoı̈de divisée par le nombre de pales. Le développement du sillage tourbillonnaire (variation du diamètre) présenté sur la figure 2.11 et figure 2.12 montre clairement pourquoi les modèles tourbillonnaires qui supposent un diamètre de sillage constant, ne permettent pas d’obtenir de bons résultats. 34 2.4 Mesures par fil chaud Les mesures PIV, permettent d’obtenir les composantes axiale et radiale de vitesse dans les plans azimutaux. Pour obtenir la composante tangentielle, les explorations sont faites à l’aide de l’anémométrie à fil chaud. L’anémométrie à fil chaud est une technique de mesure qui permet de déterminer la vitesse instantanée d’un fluide s’écoulant autour d’une sonde à fil chaud placée au sein de l’écoulement. Depuis son invention par L.V.King en 1914 [60], cet instrument de mesure est très utilisé dans le domaine de la mécanique des fluides. Il y a deux méthodes de mesure par fil chaud : l’anémomètre à température constante (CTA) [61] et l’anémomètre à courant constant (CCA) [62]. Dans nos travaux, nous avons utilisé l’anémomètre à température constante (CTA), car celui-ci a l’avantage d’avoir une plus grande bande passante ainsi qu’une fréquence de coupure plus élevée allant jusqu’à 1 MHz. 2.4.1 Principe de l’anémométrie à fil chaud Le principe de cet appareil CTA, est basé sur la loi de transfert de chaleur entre un écoulement fluide et un élément très sensible et très fin, de l’ordre de 5 µm de diamètre et de longueur 1 mm, chauffé par un courant électrique, et placé dans cet écoulement, figure (Fig.2.13). V1 E Q Pont de Wheatstone 1 V Q Servoamplificateur Fil chaud Fig. 2.13 – Principe du fil chaud. Il s’agit de maintenir à température constante un fil en tungstène traversé par un courant électrique. Si l’on met une sonde à fil chaud dans un écoulement fluide, le fil est refroidi par convection. Pour maintenir la température du fil constante, on doit augmenter 35 l’intensité électrique le traversant de façon à apporter une énergie (RΩ × I 2 ) égale à celle perdue par convection. RΩ : Résistance de fil chaud. I : Intensité du courant électrique. Le fil chaud est connecté avec un pont de Wheatstone et chauffé par un courant électrique. Une variation de la résistance du fil chaud, qui est fonction de la température, rend le pont déséquilibré. Un servoamplificateur est connecté au pont, celui-ci, maintient toujours le pont en équilibre en contrôlant le courant qui passe dans le fil chaud, par conséquence, la température du fil reste constante. Sans écoulement, on s’arrange pour avoir un courant de chauffage I. Avec un écoulement, une quantité de chaleur Q est transférée à l’écoulement et on mesure une tension E (la tension de pont de Wheatstone). Afin de maintenir la température du fil chaud constante, et donc le pont en équilibre, le courant I, renforcé par la variation la tension E, chauffe le fil chaud. Donc, une variation du courant I, et de la tension E, permet de capter les variations de la vitesse d’écoulement selon la loi de King : E 2 = G + jVlϕ E2 − G Vl = j   ϕ1 (2.1) (2.2) - ϕ ≈ 0.5 - G, j : constantes de calibrage de la pont de Wheatstone. - Vl : vitesse locale. A partir du principe du fil chaud, plusieurs catégories de fil chaud sont réalisables. Avec l’utilisation d’un fil chaud, on peu mesurer une seule composante de la vitesse. La construction de deux ou trois fils chauds ensemble, donne la possibilité de mesurer deux ou trois composantes de vitesse simultanément. La figure 2.14, représente les cas précédents. Fig. 2.14 – Trois catégories de sondes à fil chaud. 36 Un avantage d’utiliser le fil chaud est qu’il permet la mesure du taux de turbulence. De plus, les sondes à deux ou trois fils permettent la mesure des écoulements turbulents car ils ont un faible encombrement et ont un temps de réponse très rapide [63]. Ils ne perturbent presque pas l’écoulement et permettent l’analyse de fluctuations très rapides. 2.4.2 Méthodologie suivie dans l’exploration au fil chaud Le but de cette exploration est d’étudier les deux champs des vitesses, axiale et tangentielle, en amont et en aval d’une éolienne. Des explorations du sillage de la même éolienne testée précédemment par la méthode PIV, sont réalisées dans la veine d’essai, figure (Fig. 2.3), avec une vitesse axiale d’écoulement d’environ 9.4 m/s. Par contre, la vitesse de rotation de l’éolienne était approximativement de 1050 tr/min. Le fil chaud utilisé pour ce travail est de type 55 P11 (Dantec Mesurement Tchnology), il est fixé à l’aide d’un support porté par un robot d’exploration et branché avec un câble de (5 m) à la chaı̂ne de mesure. 2.4.3 Plan d’exploration Dans la veine d’essai de notre soufflerie, nous avons essayé d’explorer le plan le plus grand possible. Malgré les dimensions limitées de la veine(1.35m × 1.65m × 2m), nous avons pu explorer l’écoulement, jusqu’à 0.4D en amont du plan de rotation et 1.5D en aval environ. Le plan d’exploration a été défini de la façon suivante : on commence la mesure au niveau de l’axe de rotation, le pas selon l’axe Y (axe vertical) est toujours égal à 0.05 m, et les pas horizontaux sont variables. En amont, on a commencé à la distance 0.015 m du plan de rotation, et chaque pas est égal à 0.03 m. Par contre, en aval, on a commencé à partir de 0.015 m du plan de rotation avec un pas de 0.05 m jusqu’à la huitième mesure, puis avec un pas égal à 0.1 m pour les trois dernières mesures, figure 2.15. Pour chaque rayon de mesure choisi, le fil chaud enregistre sur chaque tour de rotation 800 valeurs pour la vitesse, ces valeurs sont synchronisées azimutalement par rapport aux pales grâce à un capteur optique qui envoie un signal de synchronisation à chaque passage d’une cible collée sur le moyeu, figure 2.16. L’acquisition est répétée des centaines de fois pour obtenir la vitesse moyennée. 37 0.1 m 0.05 m 0.03 m 0.05 m Fig. 2.15 – Plan d’exploration du sillage autour du rotor. Fig. 2.16 – Points de mesures sur des cercles complets. 38 2.4.4 Représentation du sillage A partir de la figure 2.16, une interpolation de vitesse à partir des points obtenus est réalisée afin d’avoir un champ de vitesse continu sur tout le domaine exploré comme le présente la figure 2.17. Fig. 2.17 – Représentation du sillage. La coupe longitudinale de ce volume est représentée sur la figure 2.18, sur cette figure on observe que le déficit de la vitesse axiale est très clair derrière le rotor. On observe aussi l’intersection entre cette coupe et les tubes tourbillonnaires émis par chaque pale. 39 Fig. 2.18 – Champ de vitesse axiale dans le sillage. 2.4.5 Étude de trajectoire des tourbillons marginaux Pour visualiser les phénomènes tourbillonnaires dans le sillage étudié, nous pouvons examiner des coupes cylindriques et coupes axiales. Les coupes cylindriques sont réalisées à trois différents diamètres. D’après la figure 2.19, où le diamètre est de 540 mm, les tourbillons sont visibles seulement à coté du rotor. Parce que les tourbillons, en s’éloignant du plan de rotation, s’éloignent aussi de l’axe de rotation et leurs positions deviennent dehors du diamètre 540 mm. Avec un diamètre de sillage égal à 560 mm, on trouve que des tourbillons sont visibles en s’éloignant du rotor , figure 2.20. Par contre, pour la figure 2.21, R = 290 mm, on trouve que les tourbillons lointains sont plus clairs. On constate d’après les résultats des trois figures 2.19, 2.20 et 2.21 que les tourbillons se détachent des extrémités des pales, prennent des trajectoires hélicoı̈dales par rapport à l’axe de rotation, et s’éloignent du celui-ci en fonction de leur distance au plan de rotation. 40 Fig. 2.19 – Champ de la vitesse axiale, R = 270 mm. Fig. 2.20 – Champ de la vitesse axiale, R = 280 mm. Fig. 2.21 – Champ de la vitesse axiale, R = 290 mm. 41 2 1 Fig. 2.22 – Vitesse axiale, Z/D = 0.32. Fig. 2.23 – Vitesse axiale, Z/D = 0.48. 4 3 Fig. 2.24 – Vitesse axiale, Z/D = 0.62. Fig. 2.25 – Vitesse axiale, Z/D = 0.8. 5 6 Fig. 2.26 – Vitesse axiale, Z/D = 1. Fig. 2.27 – Vitesse axiale, Z/D = 1.2. 42 8 7 Fig. 2.28 – Vitesse axiale, Z/D = 1.3. Fig. 2.29 – Vitesse axiale, Z/D = 1.4. D’autre part, nous avons présenté plusieurs sections axiales (normales à l’axe de rotation) successives du sillage en fonction de leurs distances au plan de rotation pour montrer aussi les trajectoires des tourbillons au fur et à mesure de l’éloignement du plan de rotation. A partir de la figure 2.22, à la distance Z/D = 0.32, on trouve les tourbillons qui quittent les extrémités des pales. Nous avons indiqué un tourbillon par une flèche, celui-ci vient de se détacher de l’extrémité d’une pale. Sur la figure 2.23, à la distance Z/D = 0.48, les tourbillons s’éloignent de l’axe de rotation avec des positions azimutales différentes, le tourbillon indiqué a pris un angle azimutal plus grand que dans le cas précédent. Dans la figure 2.24, à la distance Z/D = 0.62, on trouve que les tourbillons s’éloignent encore de l’axe de rotation. Par contre, la figure 2.25, à la distance Z/D = 0.8, montre bien que les tourbillons sont clairement d’un rayon plus grand que R. Ainsi, nous observons que le tourbillon indiqué a tourné autour l’axe de rotation de plus que 120°. Les figures 2.26, 2.27 et 2.28 montrent bien l’éloignement des tourbillons de l’axe de rotation de l’éolienne et les mouvements azimutaux. A la distance égale à Z/D = 1.4, les tourbillons ont complètement disparu parce qu’ils quittent le domaine de notre exploration, voir la figure 2.29. 43 2.4.6 Effet de l’éolienne sur l’écoulement à l’amont du rotor Pour présenter l’influence du rotor sur l’écoulement en amont, nous avons réalisé quatre sections successives à différentes distances par rapport au plan de rotation. Fig. 2.30 – Vitesse axiale à Z/D = 0.4 du plan de rotation en amont. Pour la première section, figure 2.30, à une distance égale à Z/D = 0.4, on observe que qu’il y a quatre régions de vitesses différentes à cause de l’influence du rotor éolien, la diminution de la vitesse axiale est plus grande à l’approche de l’axe de rotation de l’éolienne, parce que la corde de la pale augmente de 0.047 m à l’extrémité de pale jusqu’à 0.067 m au pied de pale où se trouve le moyeu qui fait une obstacle devant le l’écoulement. La vitesse la plus faible dans cette section est d’environ 3 m/s. En nous rapprochant du plan de rotation à une distance de Z/D = 0.25, figure 2.31, on observe l’apparition de l’influence des trois pales du rotor, et la surface où la vitesse est inférieure à 3 m/s est plus grande par rapport à la section précédente. A une distance du plan de rotation égale à Z/D = 0.125, figure 2.32, la présence de l’influence des pales devient plus nette, la vitesse est plus faible que dans le cas précédent, environ 1 m/s. 44 Fig. 2.31 – Vitesse axiale à Z/D = 0.25 du plan de rotation en amont. Fig. 2.32 – Vitesse axiale à Z/D = 0.125 du plan de rotation en amont. 45 Fig. 2.33 – Vitesse axiale à Z/D = 0.02 du plan de rotation en amont. A la distance de Z/D = 0.02 du plan de rotation, figure 2.33, la présence des pales et du moyeu est encore plus nette, l’influence des tourbillons détachés de l’extrémité de chaque pale est présente avec des régions de vitesse élevées (environ 11 m/s ). 2.4.7 Champs de vitesse tangentielle Nous avons utilisé l’anémométrie à fil chaud, pour avoir la possibilité de mesurer le champ de vitesse tangentielle qu’on n’a pas pu mesurer avec la PIV. Dans ce paragraphe, nous allons présenter le champ de vitesse tangentielle en amont et en aval de l’éolienne testée. Le vent arrive face à l’éolienne sans rotation, mais par l’influence de la rotation des pales, une vitesse tangentielle va apparaı̂tre. La figure 2.34, présente la vitesse tangentielle à distance de Z/D = 0.048 en amont du plan de rotation ; a partir de cette figure 2.34, nous observons que, la vitesse tangentielle prend des valeurs très faibles (très proche de zéro) sur la majorité de cette surface, mais dans les zones très proches des pales, la vitesse augmente jusqu’à 1.5 m/s. Après le passage du vent à travers les pales de rotor, la vitesse tangentielle va augmenter par rapport à la vitesse tangentielle en amont du plan de rotation, voir la figure 2.35. Cette augmentation est due aux influences des tourbillons attachés, des tourbillons margi46 Fig. 2.34 – Vitesse tangentielle à Z/D = 0.048 en amont du plan de rotation. naux et des tourbillons des pieds des pales qui forment une vitesse induite très importante [37]. Fig. 2.35 – Vitesse tangentielle à Z/D = 0.048 en aval du plan de rotation. La figure 2.35, représente la vitesse tangentielle à une distance de Z/D = 0.048 en aval du plan de rotation, nous observons que, au pied des pales, la vitesse tangentielle est élevée à cause de l’influence des tourbillons en pieds des pales. Par contre, sur cette figure 47 2.35, la vitesse tangentielle est aussi augment à l’extrémités des pales par l’influence des tourbillons marginaux. En s’éloignant en aval du plan de rotation, la vitesse tangentielle diminue car l’influence des tourbillons est diminuée. Sur la figure 2.36 à la distance de Z/D = 0.1 du plan de rotation, la vitesse tangentielle est moins forte que sur la figure 2.35 à Z/D = 0.048 en aval du plan de rotation. Fig. 2.36 – Vitesse tangentielle à Z/D = 0.1 en aval du plan de rotation. Enfin, la figure 2.37 représente un champ de vitesse tangentielle le long du plan d’exploration. Nous observons que la vitesse tangentielle est très importante sur les frontières du sillage, dans la zone d’intersection entre les tubes tourbillonnaires et le plan d’exploration. 2.5 Conclusion Dans ce travail expérimental, nous avons mis en œuvre la technique PIV et la technique d’anémomètrie à fil chaud pour l’exploration en soufflerie de l’écoulement et l’obtention de données quantitatives sur le champ de vitesse dans le sillage proche d’une éolienne à axe horizontal. Avec la technique PIV, la synchronisation des tirs lasers et des prises d’images par rapport à la position azimutale de la pale, a permis de visualiser l’écoulement dans un repère lié au rotor et de réaliser une reconstitution tridimensionnelle du champ de vitesse. Ainsi, la position des tourbillons marginaux a pu être localisée. Les résultats montrent que les tourbillons marginaux issus des extrémités des pales ne sont pas situés sur une surface 48 Fig. 2.37 – Champ de vitesse tangentielle. cylindrique comme le suppose la théorie tourbillonnaire linéaire [64]. Ils se déplacent vers l’extérieur en augmentant le diamètre du tube de courant comme le prévoit la théorie de Froude-Rankine [59]. Pour élargir le champ mesuré, les explorations PIV, ont été réalisés sur 6 fenêtres voisines présentant un certain chevauchement. Le traitement effectué après les mesures a permis de caler l’ensemble des résultats de ces fenêtres et d’obtenir le champ de vitesse jusqu’à 1.5D en aval le rotor. L’analyse des résultats montre l’interaction entre l’écoulement en aval du rotor et les vitesses induites par les tourbillons marginaux. Les mesures ont révélé la présence de structures tourbillonnaires importantes en aval du moyeu et près du pied des pales. La technique de l’anémométrie à fil chaud a permis d’explorer le champ de vitesse axiale et tangentielle en amont et en aval du rotor. Les essais par fil chaud, montrent bien le volume occupé par le sillage proche, les détachements des tourbillons marginaux de l’extrémité de chaque pale et leurs trajectoires qui ne sont pas situées sur une surface cylindrique. L’exploration de l’éolienne dans la veine de retour permet d’étudier le sillage lointain jusqu’à 6D, celui-ci est assez perturbé. 49 50 Chapitre 3 Développement d’un modèle hybride basé sur le concept du disque actif 51 3.1 Introduction : Dans ce chapitre nous présentons les éléments théoriques que nous avons utilisés pour la construction d’un modèle hybride basé sur le concept du disque actif de cette thèse. Nous commençons par la présentation de la méthode de l’élément de pale BEM (Blade Element Momentum), qui est utilisée pour déterminer les performances d’une éolienne et les forces appliquées par celle-ci sur l’écoulement. Ensuite, on présente le modèle du disque actif, où le rotor est remplacé par un disque d’épaisseur nulle et de diamètre égal à celui-ci du rotor. Dans la pratique, ce disque aura une épaisseur E. Des efforts aérodynamiques volumiques sont calculés à l’aide de la théorie de l’élément de pale et sont répartis sur un ensemble d’anneaux coaxiaux d’épaisseur radiale dr. L’écoulement global est calculé par la résolution numérique des équations de Navier-Stokes turbulentes moyennées. Cet écoulement fournit en particulier la vitesse locale incidente à chaque pale de l’éolienne et permet d’évaluer les efforts aérodynamiques appliqués sur le disque actif. Comme nous l’avons déjà fait remarquer, cette simplification est faite pour diminuer le coût du calcul en diminuant le nombre de cellules du maillage. De part sa simplicité géométrique, il suffit de mailler un disque actif et son environnement, le temps consacré à la réalisation des maillages est bien plus court. Nous validons notre modèle à l’aide de données expérimentales issue du NREL et des essais effectués dans la soufflerie de l’ENSAM Paris. Trois éoliennes de référence sont choisies ici. Les puissances en fonction de la vitesse du vent sont validées par les données des éoliennes NREL phase II [51] et NREL phase VI [65]. Le champ de vitesse calculé est comparé aux mesures expérimentales relevées sur une éolienne Rutland 503 testée dans la soufflerie à l’ENSAM Paris. 3.2 Méthode de l’élément de pale La méthode la plus utilisée pour calculer la charge appliquée par l’écoulement sur une éolienne et le comportement d’une éolienne, consiste à utiliser le bilan de quantité de mouvement appliqué à un élément de pale (Blade Element Momentum). Cette méthode est basée sur la division de l’écoulement en volumes de contrôle annulaires d’épaisseur dr, auxquels on applique le bilan de quantité de mouvement et d’énergie. Ces anneaux s’étendent de l’infini amont jusqu’à l’infini aval par rapport au rotor, figure 1.4, page 14. Les hypothèses principales de cette méthode sont : – que la vitesse induite dans le plan de rotation est égale à la moitié de la vitesse induite à l’infini aval. 52 – que l’on peut analyser l’écoulement par la division de la pale en nombre d’éléments indépendants et que la force d’un élément de pale est seule responsable de la variation de quantité de mouvement de l’air qui passe dans l’anneau balayé par cet élément. Cette hypothèse est acceptable si on suppose qu’aucune interaction radiale n’existe entre les écoulements qui passent dans des anneaux voisins balayés par les éléments de pales. Pratiquement, le facteur d’induction axial de l’écoulement est rarement uniforme, mais les résultats expérimentaux de Lock (1924) sur une hélice traversée par un écoulement prouvent que l’hypothèse de l’indépendance radiale est acceptable [37]. 3.2.1 Théorie de l’élément de pale Dans cette théorie, les forces agissant sur un élément de pale sont calculées en utilisant les caractéristiques de portance et de traı̂née d’un profil bidimensionnel. Ces caractéristiques sont obtenues en utilisant l’angle d’attaque déterminé à partir de la vitesse relative de l’écoulement dans le plan perpendiculaire à l’élément de pale (la vitesse W , figure 3.2). Ainsi, on néglige l’interaction entre les tubes d’écoulement correspondant aux éléments voisins [66]. On considère une éolienne avec B pales, leur rayon est R , leur corde est c, et l’angle de calage est β mesuré entre la ligne de portance nulle du profil et le plan de rotation. La vitesse de rotation de la pale est Ω, la vitesse du vent à l’infini amont est V1 . Ωr est la vitesse tangentielle de l’élément de pale tandis que l’on utilise les facteurs d’induction axial a et tangentiel a′ pour obtenir la vitesse relative de l’écoulement 3.1. D’après la 3.2 on peut écrire : W = q V12 (1 − a)2 + Ω2 r2 (1 + a′ )2 (3.1) On peut également tirer les relations suivantes : sin φ = Ωr(1 + a′ ) V1 (1 − a) , cos φ = W W (3.2) L’angle d’incidence : α=φ−β 53 (3.3) Ωrá V(1-a) 1 dr Ωr r r Ωr Fig. 3.1 – Un élément de pale au rayon local r et un anneau balayé par cet élément. W`i W`i ` W v i W Fn Fz W Fx V1 (1-a) V1 φ α β Ft -θ - Ωr Vθ Z Fig. 3.2 – Vitesses et forces agissant sur un élément de pale. − → La vitesse induite Wi a pour composantes : − → → Wi = −aV1~z − a′ Ωr− eθ 54 (3.4) La portance générée par un élément de pale d’épaisseur dr, normale à la direction de W est : dFz = 0.5ρW 2 Cz c dr (3.5) où Cz : Coefficient de portance. − → Et la traı̂née parallèle à W est : dFx = 0.5ρW 2 Cx c dr (3.6) où Cx est le coefficient de traı̂née. On peut reprojeter ces efforts afin d’obtenir les composantes axiales ou normales et tangentielles. La force axiale est donnée par : dFn = dFz cosφ + dFx sinφ (3.7) La force tangentielle est donnée par : dFt = dFz sinφ − dFx cosφ (3.8) En adimensionnant ces relations on peut obtenir le coefficient de force normale au plan de rotation : Cn = Cz cos φ + Cx sin φ (3.9) De même, le coefficient de force parallèle au plan de rotation, pour un élément de pale est : Ct = Cz sin φ − Cx cos φ 3.2.2 (3.10) Bilan de quantité de mouvement appliqué à un élément de pale L’étude que nous venons de faire et plus particulièrement l’équation (3.9) nous permet d’évaluer la résultante des forces aérodynamiques de B éléments de pales dans la direction axiale au rayon r : 1 (3.11) dN = ρW 2 Bc(Cx sin φ + Cz cos φ)dr 2 D’autre part, cette même résultante est égale à la variation de la quantité de mouvement axiale de l’air passant à travers un anneau de rayon r, d’où : 55 1 dN = 4πρV12 a(1 − a)rdr + ρ(2a′ Ωr)2 2πrdr (3.12) 2 où on a tenu compte de la contribution axiale de la variation de quantité de mouvement produite par le facteur d’induction tangentielle. On peut alors égaler les deux écritures de dN : 1 ρW 2 Bc(Cz cos φ + Cx sin φ)dr = 4πρ[V12 a(1 − a) + (a′ Ωr)2 ]rdr 2 On note que µ = Rr , ceci permet de simplifier l’expression ci-dessus : (3.13) W2 c (3.14) B (Cz cos φ + Cx sin φ) = 8π(a(1 − a) + (a′ λµ)2 )µ 2 V1 R Nous pouvons mener la même étude sur les couples et sur la variation de quantité de mouvement angulaire. Le couple exercé sur B éléments de pale au rayon r s’écrit : 1 dC = ρcB(Cz sin φ − Cx cos φ)W 2 rdr 2 (3.15) D’autre part, la variation de quantité de mouvement angulaire de l’air passant à travers un anneau est : dC = ρV1 (1 − a)Ωr2a′ r2πrdr = 4πρV1 (Ωr)a′ (1 − a)r2 dr (3.16) L’égalité des deux expression du couple dC permet d’écrire : 1 ρW 2 Bc(Cz sin φ − Cx cos φ)rdr = 4πρV1 (Ωr)a′ (1 − a)r2 dr (3.17) 2 Une méthode BEM classique utilise le système de deux équations (3.13) et (3.17) à deux inconnues a et a′ . Ce système, généralement résolu à l’aide d’une méthode de Newton, permet d’obtenir les caractéristiques locales de l’éolienne. A partir de ces équations, on trouve que : a= a′ = 1 2 ( 4FPσCsinn φ ) + 1 1 φ cos φ ( 4FP sin ) σCt −1 (3.18) (3.19) où la correction de Prandtl FP et la solidité de l’élément de pale σ, sont déterminables selon les équations suivantes [67, 68] : FP = B R−r 2 arccos(e(− 2 r sin φ ) ) π 56 (3.20) Bc (3.21) 2πr Il suffit ensuite d’intégrer ces efforts et couples élémentaires le long du rayon de la pale pour obtenir les performances globales de l’éolienne. Ainsi, le couple total a pour expression : σ= Cmeca = Z Rmax Rmin B B ρcCt W 2 rdr = ρc 2 2 Z Rmax Ct W 2 rdr (3.22) Rmin Ce couple permet aussi de calculer la puissance mécanique sur l’arbre du rotor, en intégrant la puissance élémentaire : B B dPmeca = dCmeca Ω = ( ρcCt W 2 rdr)Ω = ρcCt W 2 rΩdr 2 2 Pmeca = Z Rmax Rmin B B ρcCt W 2 rΩdr ⇒ Pmeca = ρcΩ 2 2 Z Rmax Ct W 2 rdr (3.23) Rmin La puissance électrique est donnée par une équation caractéristique de l’éolienne NREL phase II [51] : Pelec = 0.9036 · Pmeca − 0.847 (3.24) Nous venons ici de rappeler les principaux éléments de la méthode des éléments de pale. Nous allons maintenant l’appliquer au cas du calcul des performances du rotor de l’éolienne NREL Phase II. 3.3 Application de la méthode de l’élément de pale Nous avons appliqué la méthode de l’élément de pale sur une éolienne NREL phase II, cette éolienne a trois pales de diamètre égal à 10.1 m, le diamètre du moyeu est égal à 0.75 m, la corde des pales est constante c = 0.458 m, les pales ne sont pas vrillées, et l’angle de calage est fixé à 12 degrés. L’essai est fait à une vitesse de rotation égale à 72 tr/min, le rotor est fixé sur une nacelle, et cet ensemble est situé à une hauteur de 16.8 m [51]. 57 3.3.1 Algorithme de calcul La figure 3.3 présente l’algorithme de calcul utilisé pour calculer les forces appliquées par le rotor sur l’écoulement à l’aide de la théorie de l’élément de pale. A partir de celleci, on observe que le calcul a commencé en connaissant trois valeurs, la vitesse du vent à l’infini amont V1 , le facteur d’induction axial a et tangentiel a′ de l’écoulement. Ici on suppose que a = a′ =0, puis on prend un élément de pale situé au rayon local r, on calcule les caractéristiques du profil sur celle-ci, ensuite on détermine les nouvelles valeurs de a et a′ , donc la différence entre a au départ de notre calcul et la nouvelle valeur de a est déterminable, si cette différence est supérieure de la différence souhaitée, on recommence le calcul jusqu’à convergence. La même méthode permet d’obtenir a′ . Ensuite on applique cette opération sur n éléments de pale à partir du rayon inférieur jusqu’au rayon supérieur de la pale. Enfin, nous avons les caractéristiques propres des n éléments de pales qui nous permettent de calculer la puissance du rotor. 3.3.2 Discrétisation de la pale et critère de convergence Le calcul par la méthode de l’élément de pale impose une discrétisation du domaine de calcul. Afin de choisir une discrétisation satisfaisante, nous allons évaluer la puissance de l’éolienne obtenue en fonction de la finesse de la discrétisation choisie. On effectue cette étude sur le cas de l’éolienne NREL s809 phase II. Son rayon d’extrémité est égal à 5.05 m, le rayon du moyeu (le pied de la pale) est égal à 0.75 m, donc la longueur de la pale est égale à 4.3 m. Nous l’avons divisée en n parties, avec n = 20, 50, 100, 250, 500, 700, 800 et 1000. Nous testons ici uniquement la méthode de calcul représentée sur la 3.3 sans aucun couplage avec un logiciel externe, la vitesse incidente est imposée directement par l’utilisateur. Nous observons qu’à partir de 700 éléments de pale, dr = 0.006143 m, l’évolution de la puissance en fonction de la vitesse du vent ne varie plus. Dans la suite de ce travail, nous choisissons n = 1000 afin de s’assurer que la discrétisation spatiale est suffisante. L’algorithme itératif permettant d’obtenir les coefficients d’induction a et a′ doit également être étudié de manière à évaluer le critère de convergence requis pour arrêter l’algorithme. On quitte la boucle itérative sur a et a′ lorsque la différence de ces valeurs entre deux itérations successives |an+1 − an | et |a′ n+1 − a′ n | est inférieure à une valeur prescrite notée ǫ. On teste notre algorithme sur l’éolienne NREL s809 Phase II pour des valeurs de ǫ comprises entre 10−1 et 10−5 . Les résultats de cette étude sont représentés sur la figure 3.5. A partir de ǫ < 10−3 , le calcul de la puissance en fonction de la vitesse 58 Initialisation des V1, a et a′ Pour R min ≤ r ≤ R max Calculer les angles α et φ Lecture Cx et Cz Calculer C t et C n Calculer FP Calculer an+1et a′ n+1 n | a n+1- a |≤ ε ? n+1 n | a′ - a′ |≤ ε ? Non Oui r = R max ? Non Oui Calculer la puissance Afficher les résultats Fin Fig. 3.3 – Algorithme de calcul selon la théorie de l’élément de pale. 59 Fig. 3.4 – Influence du nombre de divisions de la pale sur l’évolution de la puissance, NREL phase II. du vent ne change plus. Le critère de convergence appliqué dans tous nos calculs est fixé à ǫ = 10−5 . Les deux paramètres évoqués ci-dessus (discrétisation spatiale et critère de convergence) ayant été validés. Nous présentons sur la figure 3.6, les résultats produits par la méthode de l’élément de pale en comparant ses résultats à la puissance électrice mesurée expérimentalement sur l’éolienne NREL s809 Phase II. L’accord entre les deux résultats est satisfaisant étant donné que nous avons ici utilisé la méthode de l’élément de pale avec un champ incident imposé par l’utilisateur. 60 Fig. 3.5 – Influence de la précision de calcul des coefficients d’induction sur l’évolution de la puissance, éolienne NREL phase II. Fig. 3.6 – Comparaison des puissances calculées par la méthode de l’élément de pale et les données expérimentales pour l’éolienne NREL phase II. 61 3.4 Modèle de disque actif En mécanique des fluides, le disque actif est défini comme une surface de discontinuité où des forces de surface agissent sur l’écoulement. Ce modèle est extrêmement simplifié et repose sur les hypothèses suivantes [59, 44] : V1 Vd P+ D V2 P- Fig. 3.7 – Modélisation du rotor éolien par un disque actif. - La géométrie du rotor est effacée et ce dernier n’est représenté que par un disque d’épaisseur nulle de diamètre D. - Le fluide est incompressible, non visqueux et non pesant. - Les vitesses V1 à l’infini amont, Vd dans le plan du disque et V2 dans la veine à l’infini aval sont uniformes et axiales. - L’énergie spécifique de l’écoulement comporte deux parties : cinétique et potentielle de pression. La vitesse axiale dans le plan de rotation est définie en fonction de la vitesse à l’infini amont par l’introduction du facteur d’induction axial a , soit : Vd = (1 − a)V1 3.4.1 (3.25) Equation de continuité L’application de l’équation de continuité permet d’écrire : ρA1 V1 = ρAd Vd = ρA2 V2 62 (3.26) 3.4.2 Bilan de quantité de mouvement Quand le vent passe dans le tube de courant comportant le disque actif, il y a un changement de vitesse égal à (V1 − V2 ), et le taux de variation de quantité de mouvement est égal à la somme des efforts extérieurs appliqués. Comme le tube de courant est complètement entouré par le vent à la pression atmosphérique, les forces à l’origine du changement de quantité de mouvement viennent uniquement de la différence de pression créée par le disque actif : (Pd+ − Pd− )Ad = (V1 − V2 )ρAd Vd Comme Vd = V1 (1 − a) Il vient alors : (Pd+ − Pd− )Ad = (V1 − V2 )ρAd V1 (1 − a) (3.27) Pour obtenir la différence de pression (Pd+ −Pd− ), nous utilisons l’équation de Bernoulli entre l’infini amont et le disque et entre le disque et l’infini aval : A l’amont, 1 1 (3.28) P1 + ρV12 = Pd+ + ρVd2 2 2 A l’aval, 1 1 P1 + ρV22 = Pd− + ρVd2 (3.29) 2 2 D’où : 1 (3.30) ∆P = Pd+ − Pd− = ρ(V12 − V22 ) 2 A l’aide des équations (3.25), (3.27) et (3.30) on trouve que : V2 = (1 − 2a)V1 (3.31) En comparant les équations (3.25) et (3.31) on trouve que la vitesse induite dans le plan du rotor est égale à la moitié de la vitesse induite à l’infini aval. 3.4.3 Coefficient de puissance D’après les équations (3.27) et (3.31), on exprime la force appliquée sur le disque actif : Ff orce = (Pd+ − Pd− )Ad = 2ρAd V12 a(1 − a) La puissance transmise au disque est : 63 (3.32) P = Ff orce Vd = 2ρAd V13 a(1 − a)2 (3.33) Le coefficient de puissance est défini par le rapport entre la puissance transmise au disque actif et une valeur de référence correspondant à la puissance du vent amont traversant une surface égale à celle du disque actif : Cp = 2ρAd V13 a(1 − a)2 1 ρV13 Ad 2 Donc : Cp = 4a(1 − a)2 3.4.4 (3.34) Limite de Betz L’équation (3.34) montre que le coefficient de puissance dépend du facteur d’induction axial a. La valeur maximale de Cp est déterminée par : dCp = 4a(1 − a)(1 − 3a) = 0 da soit a= 1 3 Ce qui correspond à, 16 = 0.593 (3.35) 27 Cette valeur est appelée la limite de Betz et montre la limite supérieure théorique de la puissance que l’on peut extraire du vent incident avec une éolienne. Cpmax = 3.5 Description du modèle hybride L’épaisseur d’un disque actif est généralement considérée nulle [33, 69]. Toutefois, l’introduction de notre modèle hybride dans un code industriel permet un meilleur contrôle des efforts appliqués localement si l’on considère des efforts répartis sur un volume et non sur une surface. En conséquence, nous utilisons ici un disque actif d’épaisseur E, comme dans la figure 3.8. Pour choisir l’épaisseur E, nous avons calculé par ce modèle hybride la puissance de l’éolienne NREL phase II en fonction de l’épaisseur E pour une vitesse à l’infini amont de 10 m/s. A partir de la figure 3.9 on observe que, l’épaisseur du modèle hybride qui donne 64 Y R Moyeu Z θ dr ∆θ r X E Fig. 3.8 – Le modèle du rotor. la puissance la plus proche de la puissance expérimentale (6.3 kW) pour l’éolienne NREL phase II est de manière assez prévisible E = c × sin(β). Pour l’éolienne en question, on a c × sin(β) = 0.208 × c = 0.09522 m. Cette épaisseur du modèle hybride est très proche de l’épaisseur du volume cylindrique réel E´balayé par les pales qui est représenté sur la figure 3.10. Dans la suite de notre travail, on choisira donc l’épaisseur du disque actif comme suit : E = c × sin β où β est l’angle de calage et c est la corde de la pale. Ainsi, on divise le disque en un nombre d’anneaux n, chacun de longueur E et de largeur radiale de dr, figure 3.8. 65 Fig. 3.9 – Évolution de la puissance en fonction de l’épaisseur du disque E/c, la vitesse à l’infini amont est V1 = 10m/s, éolienne NREL phase II. E΄ β E Fig. 3.10 – Définition de l’épaisseur du disque. 3.6 Etude de la discrétisation du disque actif Nous étudions ici l’influence de la finesse du maillage hexaédrique utilisé pour discrétiser le disque actif. A cette fin, nous allons calculer la puissance restituée par l’éolienne pour un vent à l’infini V1 = 10 m/s avec des maillages de finesse croissante. Nous présentons l’influence de ce facteur dans le cas NREL phase II. Le rayon de disque est maillé avec 1000 mailles, le pas dans la direction azimutale est égal à ∆θ = 3.6 degrés. Le calcul de la puissance de l’éolienne, avec cette densité du maillage, la vitesse à l’infini amont étant égale à 10 m/s, donne une puissance égale à 7.51597 kW (la puissance expérimentale correspondante est égale à 6.3 kW ). Ensuite, nous augmentons la densité du maillage progressivement avec la diminution du pas dans la direction azimutal. Nous avons trouvé qu’après un pas azimutal égal à ∆θ = 0.90657 degré, la puissance obtenue par ce modèle hybride ne change pas et est 66 égale à 7.38852 kW , figure 3.11. Dans notre travail, nous avons maillé le modèle hybride DA avec une densité du maillage correspondante à ∆θ = 0.90657 degré. Fig. 3.11 – Influence de la finesse du maillage sur la puissance calculée, V1 = 10 m/s, NREL phase II. 67 3.7 Interfaçage avec le logiciel Fluent Une fonction définie par l’utilisateur UDF1 , permet de paramétrer le logiciel Fluent et de le coupler avec un calcul spécial. Dans le cas de cette étude, l’intérêt va se porter sur la possibilité de calculer une distribution volumique de force différente pour chaque maille du système et non plus une distribution commune à toutes les mailles du même système. Une UDF est une fonction programmée par l’utilisateur qui peut être liée avec le solveur Fluent pour améliorer les capacités standard du code de calcul. Les UDF sont programmées dans le langage informatique C. Elles utilisent à la fois les librairies de fonctions de C et les macros prédéfinies avec le logiciel Fluent qui permettent d’accéder aux données du solveur. Les données du solveur auxquelles accède l’UDF sont : * Coordonnées du centre de la maille. * Vecteur vitesse associé à la maille. Les paramètres utilisés comme données d’entrée par le solveur sont : *Force volumique dans la direction x. *Force volumique dans la direction y. *Force volumique dans la direction z. Ces paramètres ne peuvent être calculés que pour une maille volumique ; en conséquence, il faut adapter le modèle surfacique de disque actif en un modèle volumique. Ceci fait l’objet du paragraphe suivant, qui reprend la théorie des éléments de pale. 3.8 Répartition volumique des forces appliquées sur le disque actif A partir de la théorie de l’élément de pale, on peut démontrer qu’un élément de pale reçoit une force élémentaire par l’écoulement égale à : 1 dF~ = ρcCF W 2 dr~e 2 3 ρ : Masse volumique de l’air. [kg/m ] c : Corde du profil. [m] CF : Coefficient caractéristique du profil. ~e : Vecteur unitaire définissant la direction dF~ . 1 UDF : User-Defined Function. 68 (3.36) Cette force élémentaire multipliée par le nombre de pales peut être répartie uniformément sur un anneau de rayon r et d’épaisseur radiale infiniment petite dr. Or on peut aussi ramener cette distribution à une distribution volumique de force en choisissant un élément de volume annulaire. Cet élément de volume est situé au rayon r, et a pour épaisseur radiale dr et pour longueur axiale E. E représente l’épaisseur du cylindre 3.29, page 63 qui remplace le rotor. L’élément de volume s’écrit donc comme : dv = 2πrEdr (3.37) En divisant les deux équations (3.36) et (3.37), on obtient : 1 ρcCF W 2 dr dF~ 2 = ~e dv 2πrEdr Soit l’intensité de la force égale à : BρcCF W 2 dF~ = ~e f~ = B dv 4πrE (3.38) (3.39) En projetant cette équation sur les cordonnées x, y et z, voir la figure 3.12, et à l’aide de la figure 3.2 page 54 on trouve que : La composante de l’intensité de la force appliquée par un élément de pale sur l’écoulement dans la direction x est : BρcCt W 2 sin θ (3.40) fx = − 4πrE La composante de l’intensité de la force appliquée par un élément de pale sur l’écoulement dans la direction y est : fy = − BρcCt W 2 cos θ 4πrE (3.41) θ : Angle azimutal. [deg] La composante de l’intensité de la force appliquée par un élément de pale sur l’écoulement dans la direction z est : BρcCn W 2 4πrE Le coefficient de la force normale dans le plan de rotation est : fz = Cn = Cx sin φ + Cz cos φ 69 (3.42) (3.43) Le coefficient de la force tangentielle dans le plan de rotation est : Ct = Cz sin φ − Cx cos φ (3.44) Y dr fX Ω r θ fZ fy X Fig. 3.12 – Position angulaire d’une maille et forces agissant sur elle-ci. 3.9 Algorithme de calcul L’UDF programmée pour cette étude, calcule la force volumique au centre d’une maille (point de rayon r et d’angle d’azimut θ ), appartenant au disque modélisant le rotor, puis l’affecte à la maille dans son ensemble. Cependant, puisque l’intensité de la force volumique f est appliquée uniformément sur un volume- en l’occurrence une maille- la distribution de f sur le cylindre rotor est discrétisée et la qualité de la discrétisation dépend de celle du maillage : autrement dit, dr apparaı̂t implicitement et est fonction de la finesse du maillage du cylindre rotor. Pour chaque maille appartenant au cylindre rotor, on va effectuer la boucle qui est décrite dans la figure 3.13. 70 Lecture des coordonnées du centre de la maille. Mcentre = ( x centre ,ycentre ,z centre ) Calcul du rayon du centre de la maille. r=√x²centre + y²centre tanθ = ycentre / x centre Calcul de l’azimut du centre de la maille. Vair = |V z | Lecture de la vitesse axiale du vent sur la maille. Calcul de la vitesse tangentielle. V tan = Vy cos θ - Vx sin θ Calcul de la vitesse relative. W = √ V²air+( Ω r + V tan)² Calcul de l’angle relatif φ et de l’angle d’incidence α. tanφ = Vair / ( Ω r + Vtan ) α=φ−β Calcul des coefficients de traînée et de portance du profil, C x et Cz C x = f (α, Re), Cz = f (α, Re) Cn = Cz cos φ + Cx sin φ C t = Cz sin φ − C x cos φ Calcul des coefficients des forces dans et normale au plan de rotation, Ct et C n Calcul de l’intensité des forces sur les mailles. fx= - B ρ c Ct W² 4лrE sin θ fy = - B ρ c Ct W² 4лrE cos θ fz = B ρ c Cn W² 4лrE Fig. 3.13 – UDF algorithme. Pour la simulation du sillage par le modèle hybride basé sur le concept du disque actif, cette méthode donne accès à toutes les variables permettant de calculer plusieurs grandeurs mécaniques : 71 * Effort de traı̂née de l’éolienne. * Effort tangentiel qui s’exerce sur l’éolienne. *Couple mécanique créé sur l’arbre de l’éolienne. Ces calculs sont réalisés dans une UDF qui s’exécute au terme de chaque itération de la simulation numérique. Pour calculer les efforts, on procède en multipliant l’intensité de l’effort f par le volume de la maille considérée, puis en sommant les produits pour toutes les mailles du cylindre rotor. Pour l’effort axial ou effort de traı̂née, le calcul est : N= nmailles X (fz,i V olumemailles,i ) (3.45) i=1 Avec nmailles : le nombre de mailles constituant le cylindre rotor. fz,i : la force volumique dans la direction z pour la maille i. V olumemaille,i : le volume de la maille i. L’effort tangentiel est calculé comme suit : T = nmailles X i=1 q 2 2 · V olumemaille,i ) + fy,i ( fx,i (3.46) Et le couple a pour valeur : Couplemeca = nmailles X i=1 q 2 2 ( fx,i + fy,i · ri · V olumemaille,i ) (3.47) Avec fx,i et fy,i les forces volumiques dans les directions x et y pour la maille i. La puissance mécanique se calcule en multipliant le couple mécanique par la vitesse de rotation du rotor : Pmeca = Couplemeca · Ω (3.48) La puissance électrique est donnée par une équation caractéristique de l’éolienne (NREL phase II) [51] : Pelec = 0.9036 · Pmeca − 0.847 (3.49) Avec la puissance électrique, on dispose d’une valeur que l’on va pouvoir comparer aux résultats expérimentaux du NREL sur la même éolienne, et donc discuter de la validité du modèle étudié. Par contre, la validation pour l’éolienne NREL phase VI a été faite par rapport à la puissance mécanique. 72 3.10 Application du modèle La validation de ce modèle hybride, le modèle hybride basé sur le concept du disque actif, a été faite de deux manières différentes par rapport à trois types d’éoliennes. Tout d’abord, nous avons comparé l’évolution des puissances obtenues par ce modèle hybride en fonction de la vitesse du vent à l’infini amont, avec les données expérimentales concernant les éoliennes de type NREL phase II et VI [51, 70]. Le profil de ces éoliennes est le s809, ce profil est bien étudié à l’université de TU-Delft2 qui’a déterminé ses caractéristiques aérodynamiques [71, 72, 73], nous présentons ces caractéristiques aérodynamiques dans l’annexe A. D’autre part, nous avons comparé le sillage simulé par ce modèle hybride qui représente l’éolienne Rutland 503 avec le sillage que l’on a mesuré derrière cette éolienne dans la soufflerie de l’ENSAM-Paris. Les propriétés aérodynamiques utilisées, du profil éolien (Rutland 503), ont été obtenues numériquement. La méthode de l’élément de pale (BEM) utilisée dans ce calcul, est basée sur les polaires de profil en cas bidimensionnelle et donne des résultats réalistes dans le cas d’un angle d’attaque avant le décrochage [74, 75]. Quand cette condition n’est pas respectée, l’écoulement est tri-dimensionnel. Ici, l’écoulement décroché qui tourne avec la pale, crée un écoulement dans la couche limite le long de la pale à cause des forces centrifuges, et le long de la corde à cause des forces de Coriolis. Par conséquence, le décrochage est diminué par rapport au cas bidimensionnelle, et donc le coefficient de portance est plus élevé que dans le cas bidimensionnelle, spécialement près du pied de pale [76, 77, 78, 79]. Tous les résultats sont obtenus dans le cas d’un écoulement incompressible et stationnaire. Le modèle de turbulence utilisé est k − w − SST . 3.11 Résultats et analyse Dans cette partie de notre travail, on va présenter les résultats du modèle hybride basé sur le concept du disque actif sur trois types d’éoliennes. 3.11.1 Eoliennes comparées Les tests pour les deux éoliennes NREL phase II et VI sont faits dans une grande soufflerie située à la NASA Ames Research Center Moffett Field, California. Cette soufflerie est de grande dimensions, sa section est de 24.4 m × 36.6 m [80]. 2 TU-Delft : Delft University of Technology. 73 Fig. 3.14 – Eolienne NREL phase VI dans la soufflerie de la NASA, au centre de recherche d’Ames à Moffett Field, California, USA. La section de la soufflerie est 24.4 m × 36.6 m. - L’éolienne NREL phase VI : cette éolienne a deux pales de diamètre de 10.05 m, la corde n’est pas constante, et les pales sont vrillées. La vitesse de rotation pendant l’essai est fixée à 71.63 tours/min, le rotor est fixé sur une nacelle, et cet ensemble est situé à une hauteur de 12.192 m [65]. - L’éolienne Rutland 503 : cette éolienne a trois pales de diamètre égal à 0.5 m, l’angle de calage est constant le longe de la pale est égale à 10°, la corde de la pale varie linéairement entre 0.065 m au pied de la pale jusqu’à 0.045 m à l’extrémité, le rayon du moyeu est égale à 0.065 m, la vitesse de rotation appliquée pendant tous les essais est égale à 1050 tours/min. Les essais pour cette éolienne sont faits dans la veine d’essai de la soufflerie de l’ENSAM-paris. Cette veine est semi-guidée avec une section de 1.65 m × 1.35 m, et une longueur de 2 m. 74 3.11.2 Le domaine de calcul 10 D Entrée Bloc 1 Bloc 2 Bloc 3 Bloc 4 Bloc 5 Bloc 6 Modèle hybride Sortie Tout d’abord, le modèle hybride basé sur le concept de disque actif, qui représente l’éolienne NREL phase II, est de forme cylindrique de diamètre de 10.1 m, et d’épaisseur de 0.09522 m ( c × sinβ, figure 3.29, page 63. Ce modèle hybride contient 360000 mailles, situé à 5 D de l’entrée principale et représente le rotor, figure 3.15. Le domaine de calcul est constitué de six blocs cylindriques de plus en plus grands autour du modèle hybride, ceci afin d’obtenir un maillage assez fin autour de celui-ci et de plus en plus grand en s’éloignant, avec un nombre de mailles optimale. Le bloc 1 à un diamètre de 1.5D, et une longueur de 0.25D, ce bloc englobe le modèle hybride avec un nombre de mailles égal à 7900 mailles. Le bloc 2 englobe le bloc 1, son diamètre extérieur est égal à 2D et de longueur de 0.5D, celui-ci a 11526 mailles. Le bloc 3, englobe le bloc 2, son diamètre extérieur est égal à 3D et de longueur de 1D, celui-ci a 39250 mailles. Le bloc 4, englobe le bloc 3, son diamètre extérieur est égal à 5D et de longueur de 2D, celui-ci a 217594 mailles. Le bloc 5, englobe le bloc 4, son diamètre extérieur est égal à 7D et de longueur de 7.5D, celui-ci a 232768 mailles. Le bloc 6, englobe le bloc5, son diamètre extérieur est égal à 10D et de longueur de 20, celui-ci a 314606 mailles. Donc le total des mailles du domaine de calcul est de 1180000. 20 D Fig. 3.15 – Domaine de calcul pour le modèle hybride basé sur le concept du disque actif. 75 Pour l’éolienne NREL phase VI, le maillage est très semblable au modèle précédent, l’épaisseur du volume cylindrique qui remplace le rotor est égal 0.0364 m, le nombre total des mailles dans ce cas est égal à 360000 mailles. Dans ce cas, le bloc 1 reste avec le même diamètre extérieur de 1.5D mais avec un nombre de mailles égal à 8956. Le bloc 2, reste avec les mêmes dimensions mais avec un nombre de mailles de 14890. Par contre, pour les autres blocs, 3, 4, 5 et 6 aucun changement, et le nombre total de mailles est égal à 1120000. 76 3.11.3 Comparaison des résultats numériques et expérimentaux On va commencer nos comparaisons avec l’éolienne NREL phase II puis avec l’éolienne NREL phase VI. Fig. 3.16 – Comparaison des puissances entre le modèle hybride DA et les données expérimentales, NREL phase II. Pour l’éolienne NREL phase II, les données expérimentales obtenues par le TU-Delft [71] sont comparées aux résultats du modèle hybride, voir la figure 3.16. Les résultats de puissance en fonction de la vitesse de vent du modèle hybride sont proches des données expérimentales. Le résultat de notre modèle hybride est très semblable aux résultats trouvés par YawDyn qui a utilisé la méthode BEM dans son calcul [81]. On utilise les données expérimentales de l’éolienne NREL phase VI [70], et on les compare aux résultats obtenus par le modèle hybride basé sur le concept du disque actif, figure 3.17. On remarque à partir de cette figure que, le modèle hybride étudié donne des résultats proches des données expérimentales jusqu’à la vitesse de vent environ 10 m/s, après cette valeur de la vitesse de vent, le modèle hybride donne des résultats qui s’éloignent des données expérimentales car on arrive dans la zone de décrochage. 77 Fig. 3.17 – Comparaison des puissances entre le modèle hybride DA et les données expérimentales, NREL phase VI . 3.12 Comparaison entre le sillage mesuré d’une éolienne et le sillage calculé par le modèle hybride Nous allons comparer deux champs de vitesse axiale, le première a été mesuré en aval de l’éolienne Rutland 503. Cette exploration a été faite dans la veine d’essai de la soufflerie de l’ENSAM-Paris. Le deuxième champ comparé a été calculé par le modèle hybride DA pour la même machine. La figure 3.18, représente le champ de la vitesse axiale moyenne derrière le rotor, et sur une surface de longueur allant jusqu’à Z/D=2.1 à partir de Z/D=0.25 du plan de rotation, et de hauteur allant jusqu’à Z/D=0.9 de l’axe de rotation. On a utilisé la même surface correspondante à celle-ci dans le domaine de calcul derrière le modèle hybride, figure 3.19. La différence entre les deux champs de vitesse est donnée par la figure 3.20. 78 Fig. 3.18 – Champ de vitesse axiale moyennée expérimentale mesurée dans la soufflerie de l’ENSAM-Paris. Fig. 3.19 – Champ de vitesse axiale obtenu par le modèle hybride basé sur le concept du disque actif. 79 Fig. 3.20 – Différence entre les deux champs des vitesses axiales, modèle hybride et Essai. A partir de la figure 3.20, on observe que la différence entre les deux champs de vitesse axiale est faible et que la majorité de cette erreur est assez proche de 0. En revanche, il y a des zones où la différence est assez nette, spécialement près du rotor, où cette différence est due à l’absence de la géométrie réelle des pales dans le cas du modèle hybride. 3.13 Conclusion Dans ce travail, nous avons présenté un modèle hybride qui combine un rotor simplifié avec un solveur des équations de Navier Stokes tridimensionnel. On l’a validé sur trois types d’éoliennes : éoliennes Rutland 503, NREL phase II et NREL phase VI. Le rotor est représenté par un disque de diamètre égal au diamètre du rotor et d’épaisseur liée à la corde de la pale de l’éolienne. La validation de ce modèle hybride est faite de deux façons différentes. -Pour les performances de l’éolienne en fonction de la vitesse du vent : cette validation est faite pour l’éolienne NREL phase II et VI. Pour l’éolienne NREL phase II, on trouve qu’il y a une corrélation forte entre les données expérimentales et les résultats obtenus par le modèle hybride pour toutes les valeurs des vitesses de vent utilisées. Par contre, pour l’éolienne NREL phase VI, les résultats obtenus par le modèle hybride sont assez proches 80 des données expérimentales jusqu’à V1 = 10 m/s, au delà de cette valeur, le modèle hybride donne des résultats différents des données expérimentales à cause du décrochage mal pris en compte par le modèle DA. -Pour le champ de vitesse axiale dans le sillage : la validation est faite entre le champ de vitesse axiale obtenu par le modèle hybride et le champ mesuré dans la soufflerie de l’ENSAM-Paris en aval de l’éolienne Rutland 503. Le résultat montre que les deux champs sont assez proches. La figure 3.20 montre que la majorité des différences entre les deux champs est proche de zéro. Les différences notables sont près du rotor, car la géométrie réelle de la pale est complètement absente. En conclusion générale, les résultats obtenus par le modèle hybride basé sur le concept du disque actif ont été présentés dans ce chapitre, ils sont comparables avec les données expérimentales concernant les trois cas d’éoliennes présentées dans ce travail. La validation de ce modèle, avec les restrictions qui s’y appliquent, est considérée comme acquise. 81 82 Chapitre 4 Développement d’un modèle hybride basé sur le concept de la ligne active 83 4.1 Introduction Le chapitre précédent a présenté un modèle hybride basé sur le concept du disque actif. Les forces axiales et tangentielles dans le disque actif sont distribuées radialement avec une intensité variable mais la distribution dans la direction azimutale reste uniforme et cette distribution ne permet donc pas la représentation individuelle des pales. Pour surmonter cette limitation, et obtenir un sillage plus proche de la réalité, Sørensen et Shen [42] ont développé une représentation tenant compte du nombre des pales. Dans leur modèle, la géométrie réelle des pales est remplacée dans le solver RANS1 par des forces volumiques distribuées le long de lignes actives représentant les pales. Leurs résultats ont montré un bon accord avec la courbe de puissance mesurée pour l’éolienne Nordtank tri pales 500 kW [42]. Par la suite, Mikkelsen [43] a utilisé cette approche de ligne active avec une distribution des forces dans des plans normaux à chaque ligne active selon une répartition Gaussienne. Ensuite, Ivanell [44] a utilisé la même approche que Mikkelsen et Sørensen et a montré l’influence des paramètres qui contrôlent la distribution de Gauss sur les résultats du calcul. Il a aussi montré qu’on peut améliorer les résultats avec l’utilisation de la correction de Prandtl. Le modèle hybride présenté dans ce chapitre est basé sur le concept de la ligne active, mais chaque pale est remplacée par un volume cylindrique (Cylindre Actif- CA). Les forces exercées par un élément de pale sur le fluide sont distribuées de façon uniforme à l’intérieur de l’élément de cylindre correspondant. La validation de ce modèle hybride est faite de deux façons différentes. D’abord par rapport à l’évolution de la puissance en fonction de la vitesse du vent où des comparaisons sont faites avec l’éolienne NREL phase II [51] et l’éolienne NREL phase VI [65]. Ensuite, une autre comparaison est faite entre le sillage simulé par ce modèle hybride et le sillage de l’éolienne Rutland 503 testée dans la soufflerie de l’ENSAM-Paris. 4.2 Développement d’un modèle de cylindre actif Dans ce modèle hybride, une combinaison est faite entre la méthode de l’élément de pale et les équations de Navier-Stokes. Chaque pale est remplacée par un cylindre actif CA de longueur égale à la longueur de la pale, figure 4.1, et d’un diamètre qui sera défini suite à une étude d’optimisation spécifique. L’approche de l’élément de pale sera utilisée pour calculer les efforts exercés par la pale sur le fluide à l’aide d’un sous programme intégré au sein d’un logiciel de simulation numérique. Pour définir ces efforts, les caractéristiques aérodynamiques bidimen1 RANS : Reynolds-averaged Navier-Stokes. 84 y n r dr r θ x Ω z t Fig. 4.1 – Modèle de Cylindre Actif représentant une éolienne tripale. sionnelle des profils des pales sont utilisées conjointement avec le vent incident obtenu numériquement. Les efforts appliqués par l’élément de pale seront répartis uniformément à l’intérieur du segment de cylindre correspondant à cet élément. Il s’agit donc de forces volumiques réparties sur des cellules de calcul et représentant un terme source de quantité de mouvement. 85 n dF dFz φ dFx φ Plan de rotation β φ t α W Z Fig. 4.2 – Forces exercées par l’écoulement sur un élément de pale. La figure 4.2 représente les forces aérodynamiques agissant sur un élément de pale dr situé à la distance r de l’axe de rotation. Le vecteur W , situé dans le plan du profil, représente la vitesse de l’air relative au profil après prise en compte de la vitesse induite. A partir de la figure 4.2 on trouve que les forces élémentaires de traı̂née et de portance dFx et dFz respectivement, appliquées sur l’élément de pale, sont :   dFx = 0.5 ρ W 2 c dr Cx (4.1)   dFz = 0.5 ρ W 2 c dr Cz (4.2)   dFt = dFx cosφ − dFz sinφ = 0.5 ρ W 2 c dr (Cx cosφ − Cz sinφ) (4.3) Où Cx et Cz sont les coefficients de traı̂née et de portance du profil. Pour la suite, nous définissons un référentiel t, n, r lié à la pale tournante. r représente l’axe de la pale qui tourne avec l’axe t autour l’axe fixe n. La projection des forces élémentaires dFx et dFz sur l’axe t (dans le plan de rotation), donne la force élémentaire tangentielle dFt appliquée sur l’élément de pale dr. 86 La projection des forces élémentaires dFx et dFz dans la direction normale au plan de rotation selon l’axe n, donne la force élémentaire normale dFn appliquée sur l’élément de pale dr.   dFn = dFx sinφ + dFz cosφ = 0.5 ρ W 2 c dr (Cx sinφ + Cz cosφ) (4.4)   dFt = 0.5 ρ W 2 c dr Ct (4.5)   dFn = 0.5 ρ W 2 c dr Cn (4.6) Ct = Cx cosφ − Cz sinφ (4.7) Pour dFt et dFn , on définit les coefficients Ct et Cn tels que : A partir des équations (4.3) et (4.5) on trouve que le coefficient de la force tangentielle dans le plan de rotation, est : A partir des équations (4.4) et (4.6) on trouve que le coefficient de la force normale au plan de rotation, est : Cn = Cx sinφ + Cz cosφ (4.8) Précisons que, dFt et dFn sont des forces appliquées sur le profil par l’écoulement. Par contre, les forces appliquées sur le fluide par le profil sont −dFt et −dFn . 4.3 Répartition des efforts dans un cylindre actif Dans le modèle étudié ici, on représente chaque pale par un cylindre actif, figure 4.1. Les forces créées par la pale sont distribuées comme des forces volumiques à l’intérieur d’un cylindre actif qui remplace la pale. Selon la méthode de l’élément de pale, on divise chaque pale et chaque cylindre actif en n éléments. D’après les propriétés aérodynamiques du profil de pale, on trouve qu’un élément de pale d’épaisseur radiale dr situé au rayon r reçoit une force élémentaire de la part de l’écoulement égale à : 1 ρ c CF W 2 dr~e 2 CF : Coefficient de force appliquée au profil. ~e : Vecteur unitaire définissant la direction de F~ . dF~ = 87 (4.9) Cette force sera supposée repartie uniformément à l’intérieur de l’élément dv du cylindre actif : dv = πd2 dr 4 (4.10) d : Diamètre du cylindre actif. En divisant l’équation (4.9) par l’équation (4.10), on trouve que l’intensité de force par unité du volume est : dF~ 2 ρ c CF W 2 f~ = = ~e (4.11) dv π d2 Pour faciliter la communication avec le logiciel de simulation de l’écoulement, on analyse dF~ dans le référentiel x, y, z utilisé dans ce dernier, voir la figure 4.1. Par projection de f sur x, y, z, on trouve que l’intensité de la force appliquée par l’écoulement sur un volume élémentaire dv, par rapport aux axes x, y et z est : Dans la direction x : 2 ρ c Ct W 2 sin θ (4.12) fx = π d2 Dans la direction y : 2 ρ c Ct W 2 cos θ (4.13) fy = π d2 Dans la direction z : 2 ρ c Cn W 2 (4.14) fz = − π d2 Ces forces sont appliquées par l’écoulement sur un élément de pale. 4.3.1 Choix du diamètre du cylindre actif Pour utiliser l’approche du cylindre actif proposée dans cette thèse, la première question à régler est celle concernant le choix du diamètre du cylindre actif qui est capable de bien représenter une pale réelle. Pour aider à faire ce choix, une étude préparatoire est faite. Elle consiste à comparer les résultats de simulation numérique pour le champ de vitesse créé par le cylindre actif et celui créé par la pale réelle. Pour y arriver, une étude spécifique par simulation numérique a été faite en 2D entre le profil s809 de l’éolienne NREL phase II, figure 4.3, et le profil du cylindre actif constitué par un cercle de diamètre d, figure 4.4. 88 Domaine de calcul La figure 4.3 représente le domaine de calcul autour du profil s809, celui-ci a une longueur égale à 120 c (c : la corde de la pale) et une largeur de 80 c. Le profil s809 est situé à 40 c de l’entrée, ce profil est divisé en deux parties, extrados et intrados, chaque partie est maillée avec 125 mailles resserrées au bord de fuite et au bord d’attaque. Un maillage de couche limite est plaqué sur l’intrados et l’extrados du profil, l’épaisseur de la première maille de cette couche limite est égale à 0.00065 c de manière à obtenir des valeurs de y + de l’ordre de 30. La taille des mailles augmente en s’éloignant de l’extrados et l’intrados selon une progression géométrique de raison 1.15. Les lignes AB et BF sont maillées avec 125 mailles, les lignes AJ, BC et FD sont maillées avec 250 mailles chacune. Les trois surfaces BFDC, BCJA et FOA sont maillées avec 31250 mailles chacune. Le maillage utilisé ici est structuré et quadrangulaire. Fig. 4.3 – Domaine de calcul autour du profil s809. 89 A D J B C Fig. 4.4 – Domaine de calcul autour d’un cercle équivalent au profil s809. La figure 4.4, représente le domaine de calcul autour d’un cercle destiné à remplacer l’influence du profil s809. Afin de pouvoir comparer les solutions numériques produites dans ces deux situations, on génère des maillages ayant le même nombre de mailles et la même taille caractéristique. Le domaine maillé est un cercle de diamètre de 40 c. Les lignes AB et BC sont maillées par 125 mailles. Les lignes ADC et AJC sont maillées avec 250 mailles chacune. Les deux surfaces identiques ABCD et ABCJ sont maillées avec 31250 mailles chaque avec un maillage structuré quadrangulaire. Validation de maillage Pour justifier la construction des maillages présenter précédemment sur le domaine de calcul, nous allons calculer les coefficients de traı̂née et portance du profil s809 en fonction de l’angle d’incidence, et comparer les résultats obtenus avec des autre obtenus par l’université de Dlft2 et l’université de OSU3 [82]. 2 3 Delft University of Technology. Ohio State University. 90 Tous les calculs sont réalisés à l’aide du logiciel Fluent version 6.2.18, le modèle de turbulence choisi dans tous les cas est k − w − SST , car ce modèle de turbulence permet une prévision précise des écoulements aérodynamiques avec de forts gradients opposés de pression [83], ce modèle est plus robuste que les autres modèles, même pour des applications complexes [84, 85]. Fig. 4.5 – Coefficients de portance pour le profil s809. Le nombre de Reynolds utilisé dans cette étude est égal à 106 , ce choix correspond à la valeur utilisée par le NREL, l’université de Delft et l’université d’Ohio pour calculer les coefficients de traı̂née et de portance du profil s809 [86, 52, 82]. Nous avons calculé les coefficients de portance et de traı̂née du profil s809 par un code de simulation numérique, les résultats sont présentés sur les figures (4.5 et 4.6). A partir de ces deux figures, on observe que les résultats de notre calcul sont assez proches des résultats données par OSU et par Delft pour les deux coefficients, portance et traı̂née. Nous pouvons donc valider nos paramètres de calcul. Calcul des champs de vitesse axiale autour du profil s809 et du cercle équivalent Les simulations effectuées pour trouver le cylindre actif équivalent à la pale réelle sont faites pour quatre angles de calages différents, 0, 6, 12 et 18 degrés, et plusieurs valeurs d/c du diamètre du cercle entre 0.05 et l. Dans notre calcul, nous avons calculé les forces agissant sur le profil s809, ensuite, nous avons appliqué ces forces sur le disque qui remplace le profil. Nous présentons les 91 Fig. 4.6 – Coefficients de traı̂née pour le profil s809. comparaisons sur des domaines de calcul de forme carrée de longueur égale à six cordes pour bien montrer le comportement du sillage. Fig. 4.7 – Champ de vitesse axiale autour du profil s809, l’angle de calage est β = 6°. 92 La (Fig.4.7) présente le champ de vitesse axiale autour du profil s809 avec un angle de calage égal à 6°. On observe bien l’augmentation de la vitesse axiale sur l’extrados du profil et la diminution de cette vitesse sur l’intrados. Les variations importantes du champ de vitesse se concentrent clairement dans une région proche du profil et dans son sillage. Des comparaisons ont été faites pour 13 valeurs de d/c et quatre angles d’incidence, soit 52 cas de calcul (voir les annexes B, C, D et E). Nous allons présenter trois cas, d = 0.05 c, 0.4 c et 1 c. Fig. 4.8 – Champ de vitesse axiale autour d’un cercle de diamètre de 0.05 c, β = 6°. A partir des figures 4.8, 4.9 et 4.10, on observe qualitativement les mêmes variations de vitesse autour des cercles qu’autour du profil s809. Afin de mieux quantifier les écarts entre les vitesses produites par le profil s809 et ses cercles équivalents, nous allons tracer des profils de vitesse à une corde en amont du centre de chaque cercle. C’est en effet dans cette zone que sont prises les vitesses de référence utilisées dans les différents modèles de rotors équivalents. La figure 4.11 présente la vitesse axiale extraite du calcul complet autour des cercles équivalents comparée à la vitesse axiale du profil s809, sur la figure 4.12 nous avons tracé la vitesse tangentielle. Ces résultats montrent que les profils de vitesse à une corde à l’amont du centre des cercles sont très peu sensibles au diamètre du cercle équivalent utilisé si bien que la puissance globale de l’éolienne sera la même pour tous les cas considérés. En revanche, le sillage et le champ proche sont fortement dépendants du diamètre utilisé. 93 Fig. 4.9 – Champ de vitesse axiale autour d’un cercle de diamètre de 0.4 c, β = 6°. Fig. 4.10 – Champ de vitesse axiale autour d’un cercle de diamètre de 1 c, β = 6°. Des valeurs de d = 0.3 c à d = 0.6 c donnent des résultats corrects. Le choix d’un diamètre permettant d’obtenir la même surface que celle du profil s809 nous semble donc 94 Fig. 4.11 – Comparaison entre des profils de vitesse axiale. Fig. 4.12 – Comparaison entre des profils de vitesse tangentielle. une bonne solution pour le choix du diamètre du cercle équivalent. Ce critère nous donne ici un diamètre d = 0.4 c qui sera utilisé par la suite dans toutes nos simulations. 95 4.4 Application du modèle de cylindre actif aux cas NREL s809 phase II et VI La validation du modèle hybride CA est faite de deux manières différentes, la première concerne l’évolution de la puissance en fonction de la vitesse à l’infini amont, cette validation est faite pour deux types d’éoliennes, NREL phase II et VI. La deuxième validation a été faite sur le champ de vitesse axiale de l’éolienne Rutland 503. Tout d’abord, nous allons présenter l’algorithme de calcul utilisé dans ce modèle hybride. 4.4.1 Algorithme de calcul Dans la construction de cet algorithme, les équations (4.7), (4.8), (4.12), (4.13) et (4.14) sont principalement utilisées. Pour chaque maille appartenant au cylindre actif, on effectue la boucle présentée sur la figure (Fig.4.13). Ainsi, à partir de cet algorithme, on peut calculer : - l’effort de traı̂née dans chaque cellule. - l’effort tangentiel qui s’exerce par chaque cellule. - le couple mécanique créé autour de l’axe de rotation du rotor par chaque cellule. Le calcul est fait à l’aide d’un programme écrit en langage C et intégré dans le logiciel Fluent. Ce programme calcule la force volumique au centre de chaque maille. Chaque maille dans le cylindre actif donne une force volumique différente en fonction du volume de la maille et de la distance de cette maille par rapport à l’axe de rotation r. Ces calculs sont réalisés à chaque itération de la simulation numérique. Pour calculer les efforts appliquées par une maille du cylindre actif sur l’écoulement, on multiple l’intensité de l’effort f par le volume de la maille considérée, puis la somme des forces volumiques de toutes les mailles, donne les forces globales. L’effort dans la direction axiale est : Fn = nmailles X (fz,i · V olumemaille,i ) (4.15) i=1 L’effort tangentiel est calculé comme suit : Ft = nmailles X i=1 q 2 2 ( fx,i + fy,i · V olumemaille,i ) Et le couple vaut : 96 (4.16) xcentre Lecteur des coordonnées du centre de la maille. zcentre r=√x²centre + y²centre Calcul de rayon du centre de la maille. tanθ = ycentre / x centre Calcul de l’azimut du centre de la maille. Vx , V y , Vz . Lecture de la vitesse dans le plan de réference. Calcul de la vitesse tangentielle. ycentre V tan = Vy cos θ - Vx sin θ W = √ V²z +( Ω r + V tan)² Calcul de la vitesse relative. tanφ = Vz / ( Ω r + Vtan ) α=φ−β Calcul de l’angle d’écoulement relatif φ et l’angle d’incidence α . Calcul de la traînée et de la portance du profil. Cx et Cz Cx = f (α, Re) , Cz = f (α, Re) Cn = Cz cos φ + Cx sin φ C t = Cx cos φ − C z sin φ Calcul des coefficients C t et C n Calcul l’intensité des forces volumiques dans les mailles du cylindre actif. f x= - 2 ρ c Ct W² л d2 fy = - 2 ρ c Ct W² л d2 fz = 2 ρ c Cn W² л d2 Fig. 4.13 – Algorithme de calcul pour le modèle hybride CA. 97 sin θ cos θ Couplemeca = nmailles X i=1 q 2 2 · ri · V olumemaille,i ) + fy,i ( fx,i (4.17) La puissance mécanique se calcule par la multiplication du couple mécanique par la vitesse de rotation du rotor Ω : Pmeca = Couplemeca · Ω (4.18) La puissance électrique est donnée par une équation caractéristique de l’éolienne NREL phase II [51] : Pelec = 0.9036 · Pmeca − 0.847 (4.19) Avec la puissance électrique, on dispose d’une valeur que l’on va pouvoir comparer aux résultats expérimentaux du NREL phase II et donc discuter de la validité du modèle proposé. Pour l’éolienne NREL phase VI, on compare directement la puissance mécanique. 98 4.4.2 Surface de référence Cette surface de référence nous permet de calculer les vitesses relatives et incidentes en évitant l’effet local des forces volumiques créées par le cylindre actif sur l’écoulement. Après plusieurs essais, et plusieurs visualisations des champs des vitesses, nous avons trouvé que la meilleure position pour cette surface est à une distance égale à une corde de pale en amont du plan de rotation. En raison de la périodicité du calcul (éolienne possédant trois pales), nous avons pris un tiers de cylindre correspondant à un tiers du domaine global d’étude. On divise la surface de référence en n tiers d’anneaux de largeur radiale dr, on calcule les vitesses moyennées sur chaque tiers d’anneau. La vitesse axiale moyennée calculée sur un tiers d’anneau situé au rayon r, est utilisée par la théorie de l’élément de pale pour calculer les forces agissant sur un élément de CA à partir de l’élément de pale correspondant, figure 4.14. Cylindre actif Surface de référence dr V1 c Z Fig. 4.14 – Surface de référence dans le cas d’une éolienne tripales. 99 4.4.3 Domaine de calcul Le domaine de calcul est périodique, nous avons divisé un domaine global de forme cylindrique de longueur 10D et de diamètre de 10D. Dans le cas de l’éolienne NREL phase II, on prend un tiers du domaine global de diamètre égal à 10D et de longueur 10D, figure 4.15. Le domaine est divisé en cinq parties, figure 4.16, la première de forme cylindrique de diamètre égal à la corde de la pale, l’axe de ce cylindre est le même que celui du cylindre actif (CA) qui remplace la pale, cette partie englobe le CA, le maillage de cette partie est très fin pour bien détecter le champ de vitesse autour du CA, avec un nombre de mailles égal à 39709. Fig. 4.15 – Domaine de calcul. La deuxième partie du domaine a une forme cylindrique de longueur égale à 2 cordes, parallèle à la direction de l’axe de rotation, de diamètre égal à 1.2D, la densité du maillage de cette partie est plus faible que pour la première, le nombre de mailles est de 42055. La troisième partie est de la même forme que la deuxième, mais avec une longueur de 10 cordes et un diamètre égal à 1.6D avec un nombre de mailles de 60663. La quatrième partie a la forme d’un grand cylindre de diamètre 2D, et de longueur 1.8D, le nombre de mailles de ce bloc est égal à 22078. La cinquième partie a une longueur égale à 10D et un diamètre égal à 10D, cette partie est maillée avec une faible densité de maillage, le nombre de mailles est égal à 300164. Le cylindre actif qui remplace la pale est placé à 100 Bloc 1 Bloc 2 Bloc 3 Bloc 4 Bloc 5 5D 10D Fig. 4.16 – Section longitudinale du domaine de calcul montrant les cinq blocs autour du cylindre actif. trois diamètres de l’entrée du domaine, le nombre de mailles du cylindre actif est égal à 20695 mailles, figure 4.16. Donc, le nombre total de mailles est égal à 482504 mailles. La densité du maillage diminue au fur et à mesure que l’on s’éloigne du cylindre actif. Dans le cas de l’éolienne NREL phase VI, le domaine de calcul a les mêmes dimensions que dans le cas précédent avec la même structure des blocs, 10D de longueur et 10D de diamètre et 5 blocs. L’éolienne NREL phase VI a deux pales, donc on va prendre dans le calcul périodique un demi-domaine, le diamètre du cylindre actif qui remplace la pale est égal à 0.4 fois la corde moyenne de la pale, parce que la corde n’est pas constante le long de la pale. Le maillage est analogue au cas précédent et compte 700000 cellules. 4.5 Résultats et analyse Le modèle hybride CA est introduit dans le logiciel Fluent 6.2.16. Tous les résultats sont obtenus dans le cas d’un écoulement incompressible et stationnaire, utilisant le modèle de turbulence k − w − SST . Le modèle hybride est appliqué au cas des éoliennes NREL phase II et VI. La figure 4.17, représente les données expérimentales et les résultats obtenus par le modèle hybride pour la puissance de l’éolienne NREL phase II. On trouve qu’il y a une bonne corrélation pour les vitesses de vent faibles, jusqu’à 12m/s, au-delà de cette valeur, les résultats du 101 modèle hybride s’éloignent des données expérimentales jusqu’à la valeur de 20m/s. Les résultats de notre modèle hybride est très semblable aux résultats du modèle de YawDyn [81]. Fig. 4.17 – Comparaison des puissances, modèle hybride CA, NREL phase II. Par contre, la figure 4.18 représente la comparaison entre les résultats obtenus pour le modèle hybride et les données expérimentales de l’éolienne NREL phase VI. Sur cette figure, on trouve qu’il y a une bonne corrélation des résultats jusqu’à la vitesse du vent de 9 m/s, après cette valeur, on arrive dans la zone de décrochage et les résultats sont nettement plus éloignés. 102 Fig. 4.18 – Comparaison des puissances, modèle hybride CA, NREL phase VI. 4.6 Comparaison simulation/ expérience du sillage d’une éolienne Rutland 503 modifiée Dans ce paragraphe, nous allons présenter la comparaison entre le sillage mesuré en aval d’une éolienne Rutland 503 tripale et le sillage calculé derrière un modèle hybride CA qui remplace celle-ci. La machine a été testée dans la veine d’essai de la soufflerie de l’ENSAM-Paris. Le calcul du modèle hybride a été induit dans le Fluent version 6.2.16 et le modèle de turbulence est k − w − SST . Dans les deux cas, la vitesse à l’infini amont égale à 9.4 m/s. La figure 4.19, représente le champ de vitesse axiale moyenne derrière le rotor. On représente la vitesse axiale produite par le modèle CA dans le même domaine géométrique, figure 4.20. 103 Fig. 4.19 – Champ de la vitesse axiale moyenne expérimentale en aval d’un éolien Rutland 503 modifiée. Fig. 4.20 – Champ de vitesse axiale obtenu par le modèle hybride CA. Sur la figure 4.21, nous présentons la différence entre les deux champs de vitesse axiale, mesurée et calculée, nous observons qu’il des zones où la différence est proche de zéro. Ainsi, il y a des zones où la différence arrive jusqu’à une valeur de 5 m/s à côté du rotor, où cette différence est due à l’absence de la géométrie réelle des pales dans le cas du modèle hybride. 104 Fig. 4.21 – Différence entre les deux champs des vitesses axiales, modèle hybride et Essai. 4.7 Conclusion Dans le travail présenté dans ce chapitre, un modèle hybride de cylindre actif est décrit pour évaluer les performances aérodynamiques d’une éolienne. Une combinaison entre la théorie de l’élément de pale et un code de Navier-Stokes tridimensionnelle est utilisée par ce modèle hybride. L’utilisation de la théorie d’élément de pale permet de calculer d’une part, les forces aérodynamiques appliquées sur l’écoulement par le cylindre actif qui remplace la pale. D’autre part, la résolution itérative des équations de Navier-Stokes tridimensionnelle donne le champ de vitesses en amont et en aval du rotor. Pour représenter le rotor, chaque pale est remplacée par un cylindre actif de longueur égale à celle de la pale et de diamètre égal à 0.4 fois la corde de pale. La validation de ce modèle hybride a été faite de deux façons différentes : - L’évolution de la puissance en fonction de la vitesse à l’infini amont a été comparée sur deux l’éoliennes. Pour l’éolienne NREL phase II, on observe que le modèle hybride proposé donne des résultats raisonnables par rapport aux résultats expérimentaux. Le cas de NREL phase IV, est plus difficile à représenter mais les résultats restent satisfaisant. - Concernant les sillages en aval le modèle CA et l’éolienne remplacée, une comparaison a été faite entre le champ de la vitesse axiale calculé en aval le modèle hybride et celui mesure en aval de la machine testée. La résultat est présenté sur la figure 4.21, on observe qu’il y a un bon accord entre les deux champs malgré la différence près du pied de la pale. 105 Le modèle hybride basé sur le concépt de la ligne active (CA) fournit des résultats qualitativement et quantitativement correct sur le cas de validation proposés. L’effort pour produire de tels résultats numériques est plus important que dans le cas du modèle CA : le maillage est en effet plus difficile à réaliser. 106 107 Chapitre 5 Interactions entre des éoliennes 108 5.1 Introduction Le comportement du sillage d’une éolienne est une question très importante qu’il faut prendre en compte lors de la construction d’un parc éolien. En général, dans le sillage d’une éolienne, la vitesse du vent diminue et le niveau de turbulence augmente. Dans ces conditions, si une éolienne se trouve dans le sillage d’une autre, elle recevra un vent qui a déjà cédé une bonne partie de son énergie et sa production d’énergie sera affectée. De plus, les phénomènes instationnaires existants dans le sillage sont à l’origine de charges dynamiques qui peuvent à long terme accélérer la fatigue des matériaux et réduire la durée de vie des pales et des autres composants sollicités [87, 88, 89]. Pour éviter ces problèmes, il faut prévoir une distance minimale entre les éoliennes. Cependant, la valeur financière du terrain et le besoin des infrastructures dans un parc éolien tendent à rapprocher les machines le plus possible. Alors, pour déterminer cette distance minimale entre les éoliennes et ainsi augmenter la rentabilité d’un parc éolien, il est important d’étudier et de prendre en compte les effets des éoliennes les unes sur les autres et leur contribution comme éléments de rugosité dans le parc éolien [90, 91, 92]. Les premières fermes offshore furent construites à Helgoland en Allemagne en 1989, à Blekinge en Suède en 1990 et à Vindeby au Danemark en 1991 [93]. Dans ces cas, les effets du sillage (diminution de la vitesse et augmentation de la turbulence en aval des machines) se propagent sur de plus grandes distances en aval de l’éolienne par rapport au cas d’un parc terrestre, ceci étant dû à la plus faible rugosité de la surface de la mer par rapport à celle de la surface de la terre [94]. Les paramètres les plus importants dans l’étude de l’interaction entre plusieurs machines sont : la disposition relative de chaque machine et la vitesse du vent à l’infini amont. A un degré moindre, on pourrait citer également l’intensité de la turbulence de l’écoulement incident ainsi que la stabilité atmosphérique (naturel, stable, instable)1 [95]. Dans cette partie de la thèse, une étude de l’interaction entre deux éoliennes identiques situées l’une derrière l’autre va être présentée. Nous montrerons l’évolution de la puissance de deux machines en fonction de la distance entre celles-ci et en fonction de la vitesse du vent amont. Nous avons vu que le modèle DA fournissait de bon résultats de port sur simplicité, le modèle DA est donc à privilégier. Cependant, il nous semble intéressant d’effectuer l’étude des deux modèles. Cette étude utilisera donc les deux modèles hybrides 1 Naturel : les particules d’air restent à leurs nouvelles places après la disparition des forces agissant sur les particules. Stable : chaque particule va revenir à sa position initiale après la disparition des forces agissant sur les particules. Instable : chaque particule va continuer son chemin même après la disparition des forces agissent sur les particules. 109 développés dans les chapitres précédents. Enfin, nous présenterons le comportement de plusieurs machines situées l’une en aval de l’autre. Les calculs seront effectués en supposant la présence d’éoliennes de type NREL phase II. 5.2 Simulation complète de l’interaction entre deux éoliennes à l’aide d’un modèle hybride La construction d’un parc éolien est souvent basée sur des résultats de calcul d’un sillage seul. Normalement, pour prendre en compte les effets combinés des différents sillages, l’hypothèse de superposition qui permet de cumuler l’effet des différents sillages est la plus utilisée [96, 97]. Lissaman [98] est le premier à avoir utilisé la ”superposition linéaire des perturbations” créées par le sillage des différentes machines présentes sur le parc. Cependant, cette méthode échoue à la détermination de grandes perturbations car elle surestime les déficits de vitesse, ce qui pourrait aboutir à un résultat absurde i.e. l’apparition de vitesses négatives lors de la superposition de nombreux sillages. Katic et al [99] ont utilisé le principe de superposition linéaire des zones de déficit de vitesse. Dans ce cas, l’effet cumulatif calculé du sillage de plusieurs machines s’avère être plus faible par rapport à la méthode utilisée par Lissaman. La méthode de Katic permet ainsi d’obtenir des résultats satisfaisants en comparaison avec des données expérimentales. Dans notre travail, nous adoptons une méthode plus directe que les méthodes précédentes, nous utilisons un domaine du calcul qui comporte deux modèles hybrides représentant les deux machines en même temps. Ceci donne une simulation numérique complète de l’interaction éolienne-éolienne à l’aide des modèles hybrides en évitant les problèmes de superposition. Ce travail a été fait avec le logiciel Fluent version 6.2.16 avec les modèles DA et CA. Tous les résultats sont obtenus dans le cas d’un écoulement stationnaire et incompressible avec un modèle de turbulence k − w − SST . Domaine de calcul Le domaine de calcul utilisé dans le cas du modèle hybride basé sur le concept de ligne active est le tiers d’un domaine cylindrique de longueur 20D et de diamètre 10D (D étant le diamètre de l’éolienne), figure 5.1, le rotor étudié est tripale. Chaque pale (ici 3) est modélisée par un volume cylindrique (Cylindre Actif, CA) de diamètre 0.4c, et de longueur égale à la longueur de la pale ; le maillage utilisé pour ce 110 volume compte 20695 mailles. Chaque cylindre actif est entouré par un volume cylindrique, appelé première partie du domaine, de diamètre égal à la corde et d’axe confondu avec celui du CA ; ce volume est maillé avec 39709 éléments. La deuxième partie du domaine a une forme cylindrique de longueur égale à 2 c dans la direction parallèle à l’axe de rotation de l’éolienne, et de diamètre égal à 1.2D ; le maillage de cette partie est plus grossier que la première partie du domaine avec 42055 mailles. La troisième partie du domaine est de la même forme que la deuxième mais avec une longueur de 10c, un diamètre de 1.6D et un nombre de mailles égal à 60663. La quatrième partie a une forme cylindrique de diamètre 2D, de longueur 1.8D et un nombre de mailles égal à 22078. Ces parties sont réalisées identiquement pour le deuxième modèle hybride qui remplace la deuxième machine. A la fin, une grande partie de longueur 20D et de diamètre 10D entoure toutes les parties précédentes et est maillée avec 724766 mailles. Donc le nombre total des mailles, qui prend en compte les deux machines modélisées, est égal à 973840 mailles. Bloc 1 Bloc 2 Bloc 3 Bloc 4 Bloc 5 Entrée V1 Premièr groupe Modèle hybride Sortie 20 D 5D Deuxième groupe Fig. 5.1 – Domaine de calcul contenant deux modèles hybrides basés sur le concept de la ligne active. Le modèle hybride qui remplace la première machine est situé à une distance 5D de l’entrée du domaine. Le deuxième modèle hybride remplaçant la deuxième machine placée en aval de la première est situé à différentes distances du premier modèle permettant ainsi d’observer l’influence de la distance entre les deux machines. Le domaine de calcul utilisé dans le cas du modèle hybride basé sur le concept de disque actif est un domaine cylindrique d’axe confondu avec l’axe de rotation de la machine de longueur 20D et de diamètre 10D, figure 5.2. La première éolienne est modélisée par un disque d’épaisseur non nul, ici 0.09522 m (0.09522 = c × sinβ, figure 3.29, page 63 et de diamètre égal au diamètre de l’éolienne. Ce disque est maillé avec 360000 éléments et est situé à une distance égale à 5D de l’entrée principale du domaine. La deuxième éolienne est modélisée par un disque de même géométrie et ayant le même maillage. Différentes 111 Entrée positions par rapport à la première éolienne seront prises pour étudier l’influence de la distance entre les deux machines. Chaque disque est ensuite entouré par un cylindre, désigné par bloc 1, d’épaisseur c et de diamètre 1.5D, maillé avec 45700 éléments. Le deuxième cylindre, bloc 2, qui englobe le bloc 1, a une longueur de 0.5D, un diamètre de 2D et un maillage utilisant 11526 mailles. Le bloc 3, englobant le bloc 2, a un diamètre égal à 3D, une longueur égale à 1D et 69250 mailles. Le bloc 4, entourant le bloc 3, possède un diamètre égal à 5D, une longueur de 2D et 217594 mailles. Enfin, le bloc 5, entourant le bloc 4, est un cylindre de diamètre 7D, de longueur 7.5D et utilisant 232768 mailles. Ces parties sont réalisées identiquement pour le deuxième modèle hybride qui remplace la deuxième machine. A la fin, une grande partie de longueur de 20D et diamètre 10D entoure toutes les parties précédentes et est maillée avec 485694 mailles. Donc le nombre total des mailles, en prend en compte les deux machines remplacées, est égal à 1700000 mailles environ. Bloc 1 Bloc 2 Bloc 3 Bloc 4 Bloc 5 Bloc 6 10 D Sortie V1 Modèle hybride Premièr groupe Deuxième groupe 20 D Fig. 5.2 – Domaine de calcul contenant deux modèles hybrides basés sur le concept du disque actif. 5.3 Effet de l’éloignement La distance entre deux machines est un élément très important qu’il faut prendre en compte lors de la construction un parc éolien. Généralement, cette distance est prise entre 112 5 et 9 D dans la direction dominante et entre 3 et 5 D dans la direction perpendiculaire à la direction dominante. Dans ce travail, nous présentons l’effet de l’éloignement sur plusieurs distances entre les deux machines (5 D, 6 D, 7 D, 8 D, 10 D et 15 D), dans la direction dominante, en utilisant les deux types de modèles hybrides étudiés dans cette thèse. Les éoliennes testées sont de type NREL phase II avec un profil des pales de type s809. 1- Étude de l’éloignement par le modèle hybride de Cylindre Actif (CA) : dans cette partie du travail, nous avons remplacé chaque pale du rotor par un Cylindre Actif de diamètre égale à 0.4 de la corde de la pale et de même longueur que celle-ci. Cette étude a été faite pour les distances entre éoliennes indiquées précédemment. La vitesse à l’infini amont est de 10 m/s et l’intensité de la turbulence est égale à 10 %. Les résultats obtenus sont présentés sur la figure (Fig.5.3) : 113 Distance entre les deux machines 5D 6D 7D 8D 10D 15D P1 7729.48 7730.34 7730.69 7730.61 7731.02 7730.38 P2 ((P1 − P2 )/P1 ) × 100 7304.46 5.4986 7378.62 4.5501 7452.21 3.6023 7499.77 2.9859 7599.84 1.6968 7702.43 0.3615 Fig. 5.3 – Influence de la distance entre deux éoliennes sur la puissance d’une éolienne située dans le sillage de l’autre. P1 : Puissance de la première machine.[W ] P2 : Puissance de la deuxième machine.[W ] A partir de la figure 5.3, on observe que la valeur de perte de puissance [((P1 − P2 )/P1 ) × 100] varie fortement lorsque l’espacement des machines est compris entre 5D et 10D. Au-delà d’un espacement de 10D, la variation est moins forte. On peut donc dire que l’influence d’une éolienne sur une autre située dans son sillage est fortement diminuée après d’un espacement est supérieur ou égal à 10D. La figure 5.4 montre bien cette évolution. Fig. 5.4 – Influence de la distance entre éoliennes sur la puissance, modèle hybride CA. 114 Distance entre les deux machines 5D 6D 7D 8D 10D 15D P1 7500.43 7518.11 7518.75 7519.67 7519.24 7518.76 P2 ((P1 − P2 )/P1 ) × 100 6675.86 10.99 6769.02 9.96 6845.85 8.94 6915.1 8.04 7017.44 6.673 7246.46 3.62 Fig. 5.5 – Influence de la distance sur la puissance d’une éolienne située dans le sillage de l’autre en utilisant un modèle hybride DA. 2- Étude de l’effet de l’éloignement par le modèle hybride de Disque Actif (DA) : dans ce cas, chaque éolienne est remplacée par un disque d’épaisseur liée à la corde de la pale et de diamètre égal au diamètre du rotor. Nous avons étudié l’interaction entre deux éoliennes dans les mêmes conditions que l’étude avec le modèle hybride de CA présentée précédemment. Les résultats sont sur la figure 5.5 : L’évolution de la perte de la puissance est bien représentée par la figure (Fig.5.6). Fig. 5.6 – Influence de la distance sur la puissance d’une éolienne située dans le sillage de l’autre, modèle DA. 115 Pour comparer les résultats obtenus par le modèle hybride de DA et les résultats obtenus par le modèle hybride de CA, nous avons présenté touts les résultats sur la même figure 5.7. Fig. 5.7 – Comparaison entre les deux modèles hybrides DA et CA. A partir de la figure 5.7, on trouve qu’il y a des différences importantes entre les deux modèles hybrides. Pour une distance entre les deux machines égale à 5D, la perte de puissance est égale à 11 % en utilisant le modèle hybride de DA, par contre, à la même espacement entre les deux machines, la perte de puissance est diminuée jusqu’à 5.5 % en utilisant le modèle hybride CA. Cette différence entre les deux modèles hybrides est due au fait que le sillage derrière un modèle hybride de DA est beaucoup plus long que le sillage derrière un modèle hybride de CA. Pour bien démontrer la différence entre les sillages selon le modèle hybride utilisé, nous avons présenté une comparaison entre les deux sillages obtenus (selon le modèle hybride DA et le modèle hybride CA) dans le cas d’interaction entre deux éoliennes. Le cas choisi d’interaction est celui où l’espacement entre les deux machines est égal à 10D, la vitesse à l’infini amont est égale à 10 m/s, les deux éoliennes remplacées sont du type NREL phase II. Après convergence, nous obtenons les résultats présentés dans les deux figures suivantes 5.8 et 5.9. On observe bien à partir de ces figures que, le sillage derrière le modèle hybride de CA est plus court que le sillage derrière le modèle hybrides DA, et que 116 le déficit de vitesse en aval des modèles hybride DA est plus fort que celui en aval des modèles hybrides CA. Fig. 5.8 – Champ de vitesse axiale dans un cas d’interaction, modèle hybride CA. L’explication de cette différence est que, dans le cas d’un modèle hybride de DA, les forces calculées à partir des pales, sont distribuées sur un disque de diamètre D et d’épaisseur E. Donc, dans ce cas, il y a une surface égale à [π(R2 − r2 )] qui creé donc un déficit de vitesse plus long que le CA, R est le rayon d’extrémité de la pale, r est le rayon de moyeu. Par contre, dans le cas d’un modèle hybride de CA, les forces obtenues à partir des pales ; sont distribuées dans des cylindres remplaçant les pales, dans notre cas nous avons 3 cylindres. La longueur de chaque cylindre est égal à la longueur de la pale remplacée par celui-ci, et son diamètre égal est à 0.4 c (c est la corde de la pale). Donc, la surface qui fait obstacle au vent est égale à 3(0.4c(R − r)). 117 Fig. 5.9 – Champ de vitesse axiale dans un cas d’interaction, modèle hybride DA. 5.4 Effet de la vitesse à l’infini amont Dans cette partie, nous avons étudié l’influence de la vitesse à l’infini amont. Cette étude a été faite pour deux machines situées l’une derrière l’autre avec un espacement de 10D. 5.4.1 Effet de la vitesse en amont avec un modèle hybride CA Pour montrer l’effet de la vitesse du vent en amont sur l’interaction entre deux éoliennes, nous avons choisi trois valeurs pour cette vitesse, 8 m/s, 10 m/s et 12 m/s. Nous avons réalisé ce travail sur les deux modèles présentés précédemment, le modèle de DA et le modèle de CA. La figure 5.10 représente les résultats obtenus en utilisant le modèle hybride CA pour une vitesse à l’infini amont égale à 8 m/s. On observe que la valeur perte de puissance en fonction de la distance entre les deux éoliennes est plus grande que dans le cas où la vitesse à l’infini amont est égale à 10 m/s. Par contre, la figure 5.11, représente l’effet d’une éolienne sur l’autre située dans son sillage pour plusieurs distances entre les deux avec une vitesse à l’infini amont égale à 12 118 Distance entre les deux machines 5D 6D 7D 8D 10D 15D P1 3460.99 3460.31 3460.65 3460.54 3461.36 3461.04 P2 ((P1 − P2 )/P1 )% 3130.51 9.5487 3192.84 7.7294 3254.69 5.9514 3296.94 4.7276 3378.14 2.4043 3446.66 0.4154 Fig. 5.10 – Influence de la distance sur la puissance d’une éolienne située dans le sillage de l’autre, V1 = 8 m/s, modèle hybride CA. m/s. On observe que la perte de puissance dans ce cas, est plus faible que dans les cas à 8 et 10 m/s. Distance entre les deux machines 5D 6D 7D 8D 10D 15D P1 11198.7 11211.8 11212.3 11212.9 11212.6 11212.3 P2 ((P1 − P2 )/P1 ) × 100 10944.9 2.2665 10981.1 1.9432 11036.7 1.5608 11062.7 1.3402 11124.1 0.7896 11193.7 0.1663 Fig. 5.11 – Influence de la distance sur la puissance d’une éolienne située dans le sillage d’une autre, V1 = 12 m/s, modèle hybride CA. 119 Pour comparer les résultats des trois vitesses étudiées, 8 m/s, 10 m/s et 12 m/s, nous présentons ces résultats dans la figure (Fig.5.12). Fig. 5.12 – Influence de la vitesse à l’infini amont sur l’interaction entre deux éoliennes, modèle hybride CA. La figure 5.12, récapitule tous les résultats d’interaction pour le modèle CA. On y observe la diminution de la perte de puissance avec l’éloignement ainsi que la diminution de cette même perte lors que V1 augmente. 120 5.4.2 Effet de la vitesse en amont avec un modèle hybride DA Dans cette partie, nous allons présenter l’interaction entre deux éoliennes en utilisant le modèle hybride de DA. Nous présenterons ici, l’effet de l’espacement pour deux valeurs de vitesses : 8 m/s et 12 m/s ( le cas à 10 m/s a déjà été présenté). Distance entre les deux machines 5D 6D 7D 8D 10D 15D P1 3217.91 3233.07 3233.55 3234.41 3233.91 3233.45 P2 ((P1 − P2 )/P1 ) × 100 2568.65 20.176 2646.67 18.137 32708.8 16.228 2764.52 14.52 2845.53 12.00 3016.8 6.69 Fig. 5.13 – Influence de l’espacement entre deux machines sur la puissance, V1 = 8m/s, modèle hybride DA. Distance entre les deux machines 5D 6D 7D 8D 10D 15D P1 10880.2 10878.6 10878.7 10878.7 10878.8 10878.5 P2 ((P1 − P2 )/P1 ) × 100 10276.3 5.55 10324.9 5.09 10373.9 4.64 10420.7 4.21 10494.8 3.53 10669.7 1.92 Fig. 5.14 – Influence de l’espacement entre deux machines sur la puissance, V1 = 12m/s, modèle hybride DA Les figures 5.5, Fig.5.13, Fig.5.14, représentent les résultats de l’interaction pour une vitesse à l’infini amont égale à 8 m/s, 10 m/s et 12 m/s, l’intensité de la turbulence est égale à 10%. A partir de ces figures, on trouve que la perte de puissance diminue avec l’augmentation de la vitesse à l’infini amont. La figure 5.15 représente l’évolution de la perte de puissance en fonction de la vitesse à l’infini amont et la distance entre les deux machines. 121 Fig. 5.15 – Influence de la vitesse à l’infini amont sur l’interaction entre deux éoliennes, modèle hybride DA. 5.5 Interaction entre six machines en cascade Pour finir, nous présentons les champs de vitesse axiale pour six éoliennes représentées par des modèles hybrides DA, la distance entre deux machines consécutives est égale à 5D, la vitesse à l’infini amont est égale à 10 m/s et l’éolienne est de type NREL phase II. Sur la figure 5.16 on observe que chaque machine est située dans le sillage d’une autre, et le déficit de la vitesse axiale derrière chaque machine, augmente différemment au fur et à mesure qu’on s’éloigne de l’entrée du domaine. Le déficit entre la première et la deuxième et plus grand que celui-ci entre la deuxième et la troisième et ainsi de suite. Les puissances résultantes de ces éoliennes sont sur la figure 5.17, on observe que la perte de puissance entre deux machines consécutives diminue au fur et à mesure que l’on s’éloigne de la première machine. La différence de puissance entre la première machine et la deuxième est de 10 %, la différence entre la troisième machine et la quatrième est de 5 %, la différence entre la quatrième et la cinquième est de 1.5 % et la différence entre la cinquième machine et la sixième est inférieure à 1 %. 122 Fig. 5.16 – Champ de vitesse axiale pour 6 rotors consécutifs. Fig. 5.17 – Evolution de la puissance des éoliennes situées sur le même axe de rotation, i Ri = P P. 5.6 Conclusion Dans cette partie, nous avons étudié les effets induits par la présence de plusieurs éoliennes coaxiales en cascade. Le paramètre le plus important semble être l’espacement 123 entre les éoliennes. En effet, on observe naturellement que la puissance restituée par une éolienne située dans le sillage d’une autre augmente avec la distance inter-éolienne. D’autre part, la vitesse à l’infini amont est également un facteur important. On observe que la perte de puissance d’une éolienne dans le sillage d’une autre est moins sensible lorsque la vitesse augmente. Il est à noter que nous avons également effectué des tests sur le niveau de turbulence en entrée du domaine de calcul, cette influence est très négligeable ainsi que l’ont montré Sicot et al. [100] ainsi que Chohran [101]. Ces résultats prouvent la faisabilité de l’étude proposée, nous regrettons toutefois l’absence de résultats expérimentaux accessibles pour nous fournir un point de comparaison. Une des perspectives de ce travail consisterait donc à étudier expérimentalement l’interaction entre deux éoliennes de taille réduite. 124 Chapitre 6 Conclusion générale 125 Dans ce travail, nous nous sommes intéressés aux écoulements autour des éoliennes et plus particulièrement à leur sillage, avec comme objectif le développement d’une méthode de simulation numérique capable de représenter le sillage éolien sans nécessiter un temps excessif de calcul. Ainsi, avec une telle méthode, il sera possible d’étudier les performances de plusieurs éoliennes dans un parc. Pour cela, notre travail a été orienté vers le développement de modèles hybrides couplant des modèles de rotor simplifiés, comme le disque ou le cylindre actif, à des simulations numériques de l’écoulement autour du rotor. Pour analyser les phénomènes existants dans le sillage, avoir une idée claire de la structure du sillage proche d’une éolienne, et obtenir des données pouvant servir de référence à la validation de la simulation, un travail expérimental en soufflerie a été réalisé sur une maquette d’éolienne Rutland 503 modifiée. Des explorations, synchronisées avec une pale de référence, du sillage ont été faites par la technique PIV dans différents plans azimutaux. Les mesures PIV donnent les composantes de vitesse axiale et radiale. Pour explorer le champ de vitesse tangentielle, la technique de l’anémométrie à fil chaud a été utilisée pour obtenir les champs de vitesses (tangentielle et axiale) en amont et en aval du rotor. La synchronisation des mesures, par rapport à la position azimutale de la pale de référence a permis de faire une reconstruction tridimensionnelle du champ de vitesse dans le sillage. Ces explorations ont permis de visualiser les tourbillons marginaux émis des extrémités des pales. L’intersection entre les plans d’exploration azimutaux et les tubes tourbillonnaires hélicoı̈daux a permis de visualiser les foyers tourbillonnaires et le champ de vitesse induite par chaque foyer. Ainsi, la position des tourbillons marginaux a pu être localisée et leur pas déterminé. Le diamètre du tube de courant augmente en aval du rotor à cause du ralentissement de l’écoulement créé par l’éolienne. Nos résultats montrent que les tourbillons marginaux issus des extrémités des pales ne sont pas situés sur une surface cylindrique comme le suppose la théorie tourbillonnaire linéaire ; ils se déplacent vers l’extérieur en augmentant le diamètre du tube de courant comme le prévoit la théorie de Froud-Rankine. L’examen des champs successifs dans le temps montre une fluctuation et une instationnarité de la position des noyaux des tubes tourbillonnaires dans le sillage. Cette fluctuation est due à la variation temporelle de la puissance du rotor. En effet, lors d’une augmentation de la puissance absorbée, la vitesse moyenne dans le sillage diminue conduisant ainsi à une diminution du pas des tubes tourbillonnaires. En conséquence, l’équation de continuité fait que le rayon du tube de courant augmente et déplace les tourbillons vers l’extérieur. On peut noter que la fluctuation radiale des noyaux dépend faiblement de leurs positions axiales. Par contre, leur battement axial augmente vite parce que ce battement est amplifié en fonction du nombre de pas parcourus. Ainsi, il est possible à l’aide de coupes cylindriques et des coupes normales à l’axe du sillage, d’observer 126 les tourbillons marginaux et de suivre leurs trajectoires hélicoı̈dales. Pour les travaux de simulation numérique du sillage éolien, deux modèles hybrides ont été étudiés : un modèle de disque actif (DA) et un modèle de cylindre actif (CA) couplés avec une simulation numérique basée sur la résolution des équations de Navier-Stokes. La résolution itérative des équations de Navier-Stokes donne le champ de vitesse et le calcul des forces appliquées sur le DA ou CA est faite par la méthode de l’élément de pale. La validation de ces modèles hybrides a été faite de deux manières différentes : calcul de la puissance en fonction de la vitesse du vent et calcul du champ de vitesse dans le sillage. Les puissances calculées par les modèles hybrides ont été comparées avec les données expérimentales concernant les éoliennes, NREL phase II et phase VI. Les résultats montrent qu’il y a de bonnes corrélations entre les valeurs calculées et les données expérimentales notamment par faible vent. Avec l’augmentation de la vitesse du vent, les puissances fournies par les modèles hybrides sont moins prédictives, plus particulièrement dans le cas NREL phase VI. Cette divergence est probablement due à l’entrée dans une zone de fonctionnement en décrochage. En ce qui concerne le sillage, nous avons comparé le champ de vitesse axiale mesuré en aval l’éolienne Rutland 503 avec le champ de vitesse axiale calculé en aval par les modèles hybrides. La comparaison entre les deux champs, montre que la différence est assez faible et permet de valider nos modèles. Après cette validation, une étude d’interaction entre des éoliennes a été faite à l’aide des deux modèles hybrides présentés précédemment. Les résultats montrent que le phénomène d’interaction est influencé par deux facteurs importants : l’espacement entre les machines et la vitesse du vent à l’infini amont. En effet, lorsqu’une éolienne est située dans le sillage d’une autre machine, l’augmentation de la distance entre les deux éoliennes, conduit à une augmentation de la puissance fournie par la machine à l’aval. Ce résultat s’explique facilement par le fait que le sillage de la première machine a plus de temps pour récupérer l’énergie cinétique du vent. En revanche, avec une vitesse à l’infini amont plus élevée, le phénomène de mélange entre le sillage d’une machine et du vent environnant devient plus importante. Ce mélange participe à la récupération de l’énergie cinétique par le sillage et augmente la vitesse de ce dernier. Dans ce cas, la production d’énergie de la deuxième machine est améliorée. L’étude de l’interaction de plusieurs machines en cascade avec un espacement constant entre deux machines successives, montre que la différence des puissances entre deux éoliennes successives diminue au fur à mesure que l’on s’éloigne de la première machine. La suite de ce travail peut consister à affiner les conditions de simulation pour mieux l’approcher de la réalité physique. En effet, pour développer les deux modèles hybrides étudiés dans cette thèse, on a considéré une distribution uniforme de la vitesse à l’infini 127 amont et il serait souhaitable de prendre en compte la présence du gradient de vitesse crée par la couche limite terrestre. De même, la présence du mat de l’éolienne doit être prise en compte dans la simulation. Par manque de données expérimentales sur l’interaction entre éoliennes, il serait intéressant d’étudier en soufflerie le développement du sillage et la production d’énergie en cas d’interaction entre deux éoliennes. Avec ces développements, on pourra améliorer la simulation d’un ensemble d’éoliennes placées dans parc éolien. 128 Chapitre 7 Annexe 130 Annexe A Profil s809 Le profil s809, a été développé par le NREL et a été optimisé pour améliorer la production d’énergie éolienne. Celui-ci a une base de données expérimentales bien documentée qui inclut des distributions de pression, des endroits de séparation de la couche limite, des données de traı̂née, des valeurs de traı̂née et de portance (Somers 1997, Butterfield et al 1992) [65, 73]. Fig. 1 – Coefficients de portance et de traı̂née pour l’éolienne NREL phase II, profil s809. Les cordonnées du profil s809 sont sur la figure (Fig.2), L’allrle du profil est donnée par la figure (Fig.3). Extrados x y 0.00037 0.00275 0.00575 0.01166 0.01626 0.02133 0.03158 0.03136 0.05147 0.04143 0.07568 0.05132 0.1039 0.06082 0.1358 0.06972 0.17103 0.07786 0.2092 0.08505 0.24987 0.09113 0.29259 0.09594 0.33689 0.09933 0.38223 0.10109 0.42809 0.10101 0.47384 0.09843 0.52005 0.09237 0.56801 0.08356 0.61747 0.07379 0.66718 0.06403 0.71606 0.05462 0.76314 0.04578 0.80756 0.03761 0.84854 0.03017 0.88537 0.02335 0.91763 0.01694 0.94523 0.01101 0.96799 0.00600 0.98528 0.00245 0.99623 0.00054 1 0 z 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Intrados x y 0.00140 -0.00498 0.00933 -0.01272 0.02321 -0.02162 0.04223 -0.03144 0.06579 -0.04199 0.09325 -0.05301 0.12397 -0.06408 0.15752 -0.07467 0.19362 -0.08447 0.23175 -0.09326 0.27129 -0.1006 0.31188 -0.10589 0.35328 -0.10866 0.39541 -0.10842 0.43832 -0.10484 0.48234 -0.09756 0.52837 -0.08697 0.57663 -0.07442 0.62649 -0.06112 0.67710 -0.04792 0.72752 -0.03558 0.77668 -0.02466 0.82348 -0.01559 0.86677 -0.00859 0.90545 -0.00370 0.93852 -0.00075 0.96509 0.00054 0.98446 0.00065 0.99612 0.00024 1 0 0 0 Fig. 2 – Cordonnées du profil s809. z 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Fig. 3 – Profil s809. Annexes B, C, D, E Les annexes B, C, D et E présentent les résultats de simulation du profil s809 et des cercles équivalents (d = 0.05c − 1c), qui ont permis de choisir le diamètre du cylindre actif optimal. Les forces aérodynamiques appliquées sur le profil s809 sont uniformement distribuées à l’intérieur du cercle qui remplace le profil. Calages du profil s809 étudies : 0°, 6°, 12°, 18°. Annexe B Champs de vitesse axiales autour du profil s809 et des cercles équivalentes (d = 0.05c − 1c). L’angle de calage du profil est β = 0°. Fig. 4 – Champ de vitesse axiale autour le profil s809, β = 0°. Fig. 5 – Champ de vitesse axiale autour d’un cercle de diamètre 0.05 c, β = 0°. Fig. 6 – Champ de vitesse axiale autour d’un cercle de diamètre 0.1 c, β = 0°. Fig. 7 – Champ de vitesse axiale autour d’un cercle de diamètre 0.2 c, β = 0°. Fig. 8 – Champ de vitesse axiale autour d’un cercle de diamètre 0.3 c, β = 0°. Fig. 9 – Champ de vitesse axiale autour d’un cercle de diamètre 0.4 c, β = 0°. Fig. 10 – Champ de vitesse axiale autour d’un cercle de diamètre 0.5 c, β = 0°. Fig. 11 – Champ de vitesse axiale autour d’un cercle de diamètre 0.6 c, β = 0°. Fig. 12 – Champ de vitesse axiale autour d’un cercle de diamètre 0.7 c, β = 0°. Fig. 13 – Champ de vitesse axiale autour d’un cercle de diamètre 0.8 c, β = 0°. Fig. 14 – Champ de vitesse axiale autour d’un cercle de diamètre 1 c, β = 0°. Annexe C Présentation des champs de vitesse axiales autour du profil s809 et des cercles équivalentes. L’angle de calage du profil est β = 6°. Fig. 15 – Champ de vitesse axiale autour le profil s809, β = 6°. Fig. 16 – Champ de vitesse axiale autour d’un cercle diamètre de 0.05 c, β = 6°. Fig. 17 – Champ de vitesse axiale autour d’un cercle de diamètre 0.1 c, β = 6°. Fig. 18 – Champ de vitesse axiale autour d’un cercle de diamètre 0.2 c, β = 6°. Fig. 19 – Champ de vitesse axiale autour d’un cercle de diamètre 0.3 c, β = 6°. Fig. 20 – Champ de vitesse axiale autour d’un cercle de diamètre 0.4 c, β = 6°. Fig. 21 – Champ de vitesse axiale autour d’un cercle de diamètre 0.5 c, β = 6°. Fig. 22 – Champ de vitesse axiale autour d’un cercle de diamètre 0.6 c, β = 6°. Fig. 23 – Champ de vitesse axiale autour d’un cercle de diamètre 0.7 c, β = 6°. Fig. 24 – Champ de vitesse axiale autour d’un cercle de diamètre 0.8 c, β = 6°. Fig. 25 – Champ de vitesse axiale autour d’un cercle de diamètre 1 c, β = 6°. Annexe D Présentation des champs de vitesse axiales autour du profil s809 et des cercles équivalentes. L’angle de calage du profil est β = 12°. Fig. 26 – Champ de vitesse axiale autour le profil s809, β = 12°. Fig. 27 – Champ de vitesse axiale autour d’un cercle de diamètre de 0.05 c, β = 12°. Fig. 28 – Champ de vitesse axiale autour d’un cercle de diamètre de 0.1 c, β = 12°. Fig. 29 – Champ de vitesse axiale autour d’un cercle de diamètre de 0.2 c, β = 12°. Fig. 30 – Champ de vitesse axiale autour d’un cercle de diamètre de 0.3 c, β = 12°. Fig. 31 – Champ de vitesse axiale autour d’un cercle de diamètre de 0.4 c, β = 12°. Fig. 32 – Champ de vitesse axiale autour d’un cercle de diamètre de 0.5 c, β = 12°. Fig. 33 – Champ de vitesse axiale autour d’un cercle de diamètre de 0.6 c, β = 12°. Fig. 34 – Champ de vitesse axiale autour d’un cercle de diamètre de 0.7 c, β = 12°. Fig. 35 – Champ de vitesse axiale autour d’un cercle de diamètre de 0.8 c, β = 12°. Fig. 36 – Champ de vitesse axiale autour d’un cercle de diamètre de 1 c, β = 12°. Annexe E Présentation des champs de vitesse axiales autour du profil s809 et des cercles équivalentes. L’angle de calage du profil est β = 18°. Fig. 37 – Champ de vitesse axiale autour le profil s809, β = 18°. Fig. 38 – Champ de vitesse axiale autour d’un cercle de diamètre 0.05 c, β = 18°. Fig. 39 – Champ de vitesse axiale autour d’un cercle de diamètre 0.1 c, β = 18°. Fig. 40 – Champ de vitesse axiale autour d’un cercle de diamètre 0.2 c, β = 18°. Fig. 41 – Champ de vitesse axiale autour d’un cercle de diamètre 0.3 c, β = 18°. Fig. 42 – Champ de vitesse axiale autour d’un cercle de diamètre 0.4 c, β = 18°. Fig. 43 – Champ de vitesse axiale autour d’un cercle de diamètre 0.5 c, β = 18°. Fig. 44 – Champ de vitesse axiale autour d’un cercle de diamètre 0.6 c, β = 18°. Fig. 45 – Champ de vitesse axiale autour d’un cercle de diamètre 0.7 c, β = 18°. Fig. 46 – Champ de vitesse axiale autour d’un cercle de diamètre 0.8 c, β = 18°. Fig. 47 – Champ de vitesse axiale autour d’un cercle de diamètre 1 c, β = 18°. Bibliographie [1] D. Gouriéres. Energie éolienne. Université de Dakar, 1982. [2] R. Ouziaux, P. P. Jean, and J. Perrier. Mécanique des fluides appliquée. Dunod, 1966. [3] Dodge. D. M. Illustrated history of wind power development. Darrell Dodge and TelosNet Web Development, Septembre, 2001. [4] B. Sørensen. History of, and Recent Progress in, Wind-Energy Utilization. Annual Review of Energy and the Environment, 20(1) :387–424, 1995. [5] A. Mirecki. Étude comparative de chaı̂nes de conversion d’énergie dédiées à une éolienne de petite puissance. Thèse de doctorat à l’institut polytechnique de Toulouse, Avril, 2005. [6] Joris PEETERS. Simulation of dynamic drive train loads in a wind turbine. Thèse de doctorat de l’Université de Katholieke, Belgium, 2006. [7] J. Vestergaard, L. Brandstrup, and III R. D. Goddard. A Brief History of the Wind Turbine Industries in Denmark and the United States. Energy Policy, pp. 322-327, November 2004. [8] E. Golding. The Generation of Electricity by Wind Power. E.et F. N. Spon, London, 1977. [9] A. Reeves. Wind Energy for Electric Power : A REPP Issue Brief. Renewable Energy, July, 2003. [10] V. Lauber. REFIT and RPS : Options for a harmonised Community framework. Energy Policy, 32(12) :1405–1414, 2004. [11] N. I. Meyer. Danish wind power development. Energy for Sustainable Development, 1 :18–25, 1995. [12] M. Eric Buchet. L’énergie eolienne. 5e Journée de l’Environnement et du Développement Durable - DESS Gestion de la Planète / Association GAÏA - Nice, 28 septembre 2002. 151 [13] A. R. Henderson, C. Morgan, B. Smith, H. C. Sørensen, R. J. Barthelmie, and B. Boesmans. Offshore wind energy in europe a review of the state of the art. Wind Energ. 6 :35–52 (DOI : 10.1002/we.82), 2003. [14] J. Berger. Charging Ahead : The business of Renewable Energy and What it Means for America. University of California Press, Berkley, CA, Canyon, 1997. [15] R. Harrison, E. Hau, and H. Snel. Large Wind Turbines : Design and Economics. Wiley, 2000. [16] A. Benlemilim. Promotion de l’electricité Produite à Partir des Eoliennes et Preparation d’un Cadre Reglementaire Approprie. FIER 2002, Tétouan-Maroc, 2002. [17] G. Ste-Marie. Le Développement de la Filière Eolienne au Québec et Ses Coûts. Hydro-Québec, Sections 957, 1500, 2000 et 4250 du SCFP, Juin 2005. [18] B. Chabot and L. Buquet. Le Développement de l’énergie Eolienne en France en 2005. DEWI Magazin Nr. 29, August 2006. [19] F. Jésus. Les mécanismes de projet du protocole de Kyoto . Ministère de l’économie des finances et de l’industrie, Paris, 04 juillet 2005. [20] B. Mandat. Protocole de Kyoto à la Convention-Cadre des Nations Unies sur les Changements Climatiques . FCCC/INFORMAL/83, Nations Unies 1998. [21] H. G. Beyer, B. Lange, and H. P. Waldl. Modelling Tools for Wind Farm Upgrading. Proc. European Union Wind Energy Conf, Sweden, pages 1069–1072, 20-24 may, 1996. [22] B. Lange, H. P. Waldl, A. G. Guerrero, D. Heinemann, and R. J. Barthelmie. Modelling of Offshore Wind Turbine Wakes with the Wind Farm Program FLaP. Wind Energy, 6(1) :87–104, 2003. [23] Davide Medici. Wind turbine wakes - control and vortex shedding. Technical report, KTH Mechanics, Royal Institute of Technology, Stockholm, Sweden, 2004. [24] M. Magnusson, K. G. Rados, and S. G. Voutsinas. A Study of the Flow Downstream of a Wind Turbine Using Measurments and Simulations. Wind Engineering Vol.20 N.6, 1996. [25] S. G. Voutsinas, K. G. Rados, and A. Zervos. On the analysis of wake effects in wind parks. Wind Engng, 14 :204±219, 1990. [26] K. Badreddinne, H. Ali, and A. David. Optimum project for horizontal axis wind turbines ‘OPHWT’. Renewable Energy, 30(13) :2019–2043, 2005. [27] C. Leclerc. Simulation numérique de l’écoulement tridimensionnel turbulent dans un parc éolien. National Library of Canada= Bibliothque nationale du Canada, 2000. [28] H. Mattsson and E. Djerf. Evaluation of the software program windfarm and comparisons with mesure data from Alsvik. The Aeronautical Research Institute of Sweden, FFA TN 30, 2000. [29] J. Tangler and G. Bir. Evaluation of RCAS inflow models for wind turbine analysis. NREL/TP-500-35109, February 2004. [30] F. Le Chuiton. Actuator disc modelling for helicopter rotors. Aerospace Science and Technology, 8(4) :285–297, 2004. [31] D. M. O’Brien and M. J. Smith. Analysis of rotor-fuselage interactions using various rotor models. 43rd AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, January 10-13, 2005. [32] H. Snel. Review of the present status of rotor aerodynamics. Wind Energy, 1(s 1) :46–69, 1998. [33] D. J. Sharpe. A general momentum theory applied to an energy-extracting actuator disc. Wind Energy, 7(3) :177–188, 2004. [34] R. E. Froude. On the part played in propulsion by differences of fluid pressure. Transactions of the Institute of Naval Architects, 30 :390–405, 1889. [35] W. J. M. Rankine. On the Mechanical Principles of the Action of Propellers. Institution of Naval Architects, Sweden, 1865. [36] F. W. Lanchester. A Contribution to the Theory of Propulsion and the Screw Propeller’. Transactions Institution of Naval Architects, pages 57–98, 1915. [37] T. Burton, D. Sharpe, N. Jenkins, and E. Bossanyi. Wind Energy Handbook. John Wiley & Sons, 2001. [38] J. Martinez, L. Bernabini, O. Probst, and C. Rodriguez. An Improved BEM Model for the Power Curve Prediction of Stall-regulated Wind Turbines. Wind EnergyChichester, 8(4) :385, 2005. [39] A. Ahlström. Aeroelastic simulation of wind turbine dynamics. Thèse de doctorat de Department of Mechanics, Royal Institute of Technology, SE-100 44 Stockholm, Sweden, April, 2005. [40] C. Alinot and C. Masson. Aerodynamics of wind turbines in thermally stratified turbulent atmospheric boundary layer. Proceedings of the 10th Annual Conference of the CFDSC, Ottawa, Canada, pages 553–559, 9-11 Juine 2002. [41] W. Z. Shen, Sørensen J. N., and R. Mikkelsen. Tip Loss Correction for Actuator/Navier–Stokes Computations. Journal of Solar Energy Engineering, 127 :209, 2005. [42] J. N. Sørensen and W. Z. Shen. Numerical modeling of wind turbine wakes. Journal of Fluids Engineering(Transactions of the ASME), 124(2) :393–399, 2002. [43] R. Mikkelsen. Actuator disc methods applied to wind turbines. Thèse de doctorat à la Technical University de Denmark, 2003. [44] S. S. A. Ivanell. Numerical computations of wind turbine wakes. Royal Institute of Technology, Stockholm, Sweden, April 2005. [45] F. Massouh, I. Dobrev, and M. Rapin. Numerical Simulation of Wind Turbine Performance Using a Hybrid Model. 44th AIAA Aerospace Sciences Meeting and Exhibit, pages 1–10, 2006. [46] C. Masson, A. Sma, and C. Leclerc. Aerodynamic Analysis of HAWTs Operating in Unsteady Conditions. Wind Energ, 2001. [47] C. Masson and A. Smaı̈li. Numerical Study of Turbulent Flow around a Wind Turbine Nacelle. Wind Energy-Chichester, 9(3) :281, 2006. [48] L. Sons. Re-examining the Precepts of the Blade Element Momentum Theory for Coning Rotors. Wind Energ. 9 :457–478, June 2006. [49] D. Simms, S. Schreck, M. Hand, and L. J. Fingersh. NREL Unsteady Aerodynamics Experiment in the NASA-Ames Wind Tunnel : A Comparison of Predictions to Measurements. National Renewable Energy Laboratory, USA, June 2001. [50] N. N. Sørensen and J. Michelsen. Navier-Stokes predictions of the NREL Phase VI rotor in the NASA Ames 80-by-120 wind tunnel. ASME Wind Energy Symposium, pages 94–105, 2002. [51] E. P. N. Duque and W. Johnson. Numerical Predictions Of Wind Turbine Power And Aerodynamic Loads For The NREL phase II Combined Experiment Rotor. AIAA Paper, 38, 2000. [52] D. A. Simms. Unsteady aerodynamics experiment phases II-VI test configurations and available data campaigns. National Renewable Energy Laboratory, USA, July, 1999. [53] M. Magnusson, K. G. Rados, and S. G. Voutsinas. Study of the flow downstream of a wind turbine using measurements and simulations. Wind Engineering, 20(6) :389– 403, 1996. [54] S. G. Voutsinas, K. G. Rados, and A. Zervos. Wake effects in wind parks. A new modelling approach. Proc. European Community Wind Energy Conf, LübeckTravemünde, Germany, page 444±447, 8—12 March, 1993. [55] V. P. A. R. I. DE and P. I. V OU. “Vélocimétrie et Granulométrie Laser en Mécanique des Fluides”. Institut von Karman de Dynamique des Fluides, Belgique, 2005. [56] F. Massouh, I. Dobrev1, and M. Jourieh. Etude par PIV du sillage proche d’une éolienne. Actes du Colloque FLUVISU11, 7-9 juin, ECL, Ecully, FRANCE, 2005. [57] R. Theunissen, P. Corieri, and M. L. Riethmuller. Application de la technique PTV à la modélisation d’écoulements d’aérosols dans les voies aériennes pulmonaires. 9e Congrès Francophone de Vélocimétrie Laser, Bruxelles, 14-17 septembre, 2004. [58] Most J.M. D. Honoré Susset, A. Développement et validation d’un algorithme de super-résolution par corrélation directe pour la PIV : méthodes d’accélération du traitement. 8e Congrès Francophone de Vélocimétrie Laser, Orsy, France, 17-20 septembre, 2002. [59] D. Dimitrov and V. Lazarov. Source d’energie renouvlables. université technique de sofia, 1999. [60] L. V. King. On the Convection of Heat from Small Cylinders in a Stream of Fluid : Determination of the Convection Constants of Small Platinum Wires with Applications to Hot-Wire Anemometry. Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character, 214 :373–432, 1914. [61] A. M. Savill. Fundamentals of Hot Wire Anemometry . Journal of Fluid Mechanics Digital Archive, 237 :692–693, 2006. [62] M. A. Kegerise and E. F. Spina. A comparative study of constant-voltage and constant-temperature hot-wire anemometers. Experiments in Fluids, 29(2) :154– 164, 2000. [63] S. Blidi and H. Miton. Use of the Hot Wire Anemometry for Velocity and Temperature Measurements in a Turbomachine. Journal de Physique III, 5(10) :1513–1535, 1995. [64] F. Nabah, A. Ettaouil, and N. Guennoun. Calcul du sillage d’un rotor éolien à axe horizontal calculation of wake of an horizontal axis wind rotor. Forum International sur les Énergies Renouvelables (FIER), Tétouan-Maroc, 2002. [65] M. Hand, D. Simms, L. J. Fingersh, D. Jager, J. Cotrell, M. Robinson, S. Schreck, and S. Larwood. Unsteady aerodynamics experiment phase VI : wind tunnel test configurations and available data campaigns. DRAFT to be published as and NREL Technical Report, 2001. [66] J. M. Jonkman. Modeling of the UAE Wind Turbine for Refinement of FAST-AD. Technical report, NREL National Reneweable Energy Laboratory, NREL/TP-50034755, 2003. [67] P. J. Moriarty. ACH, AeroDyn Theory Manual. National Renewable Energy Laboratory, 2005. Technical report, NREL/TP-500-36881. [68] J. N. Sørensen and C. W. Kock. Model for unsteady rotor aerodynamics. Journal of Wind Engineering and Industrial Aerodynamics, 58(3) :259–275, 1995. [69] J. N. Sørensen, W. Z. Shen, and X. Munduate. Analysis of wake states by a full-field actuator disc model. Wind Energy, 1(2) :73–88, 1998. [70] C. Lindenburg. Modelling of Rotational Augmentation Based on Engineering Considrations and Measurments. European Wind Energy Conference, London, 22-25 November, 2004. [71] D. M. Somers. Design and experimental results for the s809 airfoil. NREL/SR–4406918, National Renewable Energy Lab, Golden, CO (USA), 1997. [72] J. L. van Ingen, L. M. M. Boermans, and J. J. H. Blom. Low-speed airfoil section research at Delft University of Technology. International Council of the Aeronautical Sciences, Congress, 12th, Munich, West Germany, Proceedings, pages 401–416, Oct 1980. [73] W. P. Wolfe and S. S. Ochs. CFD Calculations of s809 Aerodynamic Characteristics. 16. American Society of Mechanical Engineers wind energy symposium, Reno, NV (United States), 6-9 Jan, 1997. [74] J. Tangler and J. D. Kocurek. Wind turbine post-stall airfoil performance characteristics guidelines for blade-element momentum methods. 43rd AIAA Aerospace Sciences Meeting and Exhibit, Nivada, Reno, October 2005. [75] R. Van Rooij, A. Bruining, and J. G. Schepers. Validation of some rotor stall models by analysis of the IEA Annex XVIII field data. Proceedings of EWEC European Wind Energy Conference Madrid, June, 2003. [76] Dr. P. K. C. et al. Investigating Three-Dimensional and Rotational Effects on Wind Turbine Blades by Means of a Quasi-3D Navier-Stokes Solver. Journal of Fluids Engineering, 122 :330, 2000. [77] P. K. Chaviaropoulos, I. G. Nikolaou, K. A. Aggelis, N. N. Sørensen, J. Johansen, M. O. L. Hansen, M. Gaunaa, T. Hambraus, et al. Viscous and Aeroelastic Effects on Wind Turbine Blades. The VISCEL project. Part I : 3D Navier-Stokes Rotor simulations. Wind Energy, 6(4) :365–385, 2003. [78] Z. Du and M. S. Selig. The effect of rotation on the boundary layer of a wind turbine blade. Renewable Energy, 20(2) :167–181, 2000. [79] H. Dumitrescu and V. Cardos. Rotational effects on the boundary-layer flow in wind turbines. AIAA Journal, 42(2) :408–411, 2004. [80] D. Simms, S. Schreck, M. Hand, L. Fingersh, J. Cotrell, K. Pierce, and M. Robinson. Plans for Testing the NREL Unsteady Aerodynamics Experiment 10-m Diameter HAWT in the NASA Ames Wind Tunnel. National Renewable Energy Laboratory, 2000. [81] G. Xu and L. N. Sankar. Computational study of horizontal axis wind turbines. Journal of Solar Energy Engineering, 122 :35, 2000. [82] JL Tangler. The Nebulous Art of Using Wind-Tunnel Airfoil Data for Predicting Rotor Performance. 21st ASME Wind Energy, Reno, Nevada, 5(2/3), 14-17 January 2002. [83] F. R. Menter, M. Kuntz R., and Langtry. Ten Years of Industrial Experience with the SST Turbulence Model, Turbulence. Heat and Mass Transfer, 4 :625–632, 2003. [84] A. G. Kanaris, Mouza K. A., and Paras S. V. Designing Novel Compact Heat Exchagers For Improved Efficiency Using a CFD Code. 1st International Conference “From Scientific Computing to Computational Engineering”, 1st IC-SCCE Athens, 8-10 September, 2004. [85] L.Davidson. An Introduction to Turbulence Models. Department of Thermo and Fluid Dynamics,Chalmers University of Technology, Götemberg, Sweden, 2001. [86] RF Ramsay, MJ Hoffman, and GM Gregorek. Effects of grit roughness and pitch oscillations on the S809 airfoil. Technical report, NREL/TP–442-7817, National Renewable Energy Lab., Golden, CO (USA), 1995. [87] A. Crespo, J. Hernandez, and S. Frandsen. Survey of modelling methods for wind turbine wakes and wind farms. Wind Energy, 2(1) :1–24, 1999. [88] A. E. Feijoo and J. Cidras. Modeling of wind farms in the load flow analysis. Power Systems, IEEE Transactions on, 15(1) :110–115, 2000. [89] N. D. Kelley and H. J. Sutherland. Damage Estimates from Long-term Structural Analysis of a Wind Turbine in a US Wind Farm Environment. Wind Energy, 1997. [90] C. Crafoord. An estimate of the interaction of a limited array of windmills. N-7713539 ; DM-16, Stockholm University.(Sweden). Meteorologiska Institutionen, 1975. [91] C. Crafoord. Interaction in limited arrays of windmills. Report DM-26, Department of Meteorology, University of Stockholm, 1979. [92] R. J. Templin. An estimation of the interaction of windmills in widespread arrays. National Aeronautical Establisshment, Laboratory Report LTR-LA-171, Ottawa, Décembre, 1974. [93] D. P. Rodmell and M. L. Johnson. The Development of Marine Based Wind Energy Generation and Inshore Fisheries in UK Waters : Are They Compatible ? Scarborough Centre for Coastal Studies, Department of Biology, University of Hull at Scarborough, Filey Road, Scarborough, YO11 3AZ, UK, 2004. [94] W. Bergström, Schlez, K. Rados, B. Lange, P. Vølund, S. Neckelmann, S. Mogensen, and L. Folkerts. ENDOW (Efficient Development of Offshore Wind Farms) : Modelling Wake and Boundary Layer Interactions. Wind Energy, 7(3) :225–245, 2004. [95] K. Rados, G. Larsen, R. Barthelmie, W. Schlez, U. Hassan, B. Lange, I. Waldl, G. Schepers, T. Hegberg, and M. Magnusson. A Comparison of Wake Model Performances in an Offshore Environment. Ewec-Conference, pages 781–784, 2001. [96] S. Frandsen and M. L. Thøgersen. Integrated fatigue loading for wind turbines in wind farms by combining ambient turbulence and wakes. Wind Engineering, 23(6) :327–340, 1999. [97] M. Magnusson and A. S. Smedman. Air flow behind wind turbines. Journal of Wind Engineering & Industrial Aerodynamics, 80(1-2) :169–189, 1999. [98] P. B. S. Lissaman. Energy effectiveness of arbitrary arrays of wind turbines. American Institute of Aeronautics and Astronautics, Aerospace Sciences Meeting, 17th, New Orleans, 15-17 Jan, 1979. [99] I. Katic, J. Højstrup, and N. O. Jensen. A simple model for cluster efficiency. Proceedings of EWEC, 86 :407–410, Rome, Italy, 1986. [100] C. Sicot, P. Devinant, T. Laverne, S. Loyer, and J. Hureau. Experimental study of the effect of turbulence on horizontal axis wind turbine aerodynamics. Wind Energ. 9 :361–370, 2006. [101] C. Cochran. The Influence of Atmospheric Turbulence on the Kinetic Energy Available During Small Wind Turbine Power Performance Testing . IEA Expert Meeting on : Power Performance of Small Wind Turbines Not Connected to The Grid, Spain, April 2002.