2014GRENI061
2014GRENI061
2014GRENI061
Présentée par
Guillaume MERCIER
Modélisation de parcs
d’hydroliennes à flux transverse
avec une méthode d’équivalence
M. Jean-Luc ACHARD
Directeur de Recherche, CNRS, Président
M. Jacques-André ASTOLFI
Maitre de conférence, Ecole Navale Brest, Rapporteur
M. Michel VISONNEAU
Directeur de Recherche, CNRS, Rapporteur
M. Jean-Luc HARION
Professeur Ecole des mines de Douai, Examinateur
M. Christian PELLONE
Chargé de Recherche, CNRS, Directeur de thèse
M. Thierry MAITRE
Maitre de Conférence, Grenoble INP, Co-Directeur de thèse
M. Christophe PEYRARD
Ingénieur de recherche, Centre de R&D EDF, Invité
M. Vincent GUENARD
Ingénieur de recherche, ADEME, Invité
2
Remerciements
3
4
à Thierry pour les nombreux débats et les discussions sur ces fameuses hydroliennes, mais
également pour m’avoir permis d’avoir une belle expérience de l’enseignement, passionnante
et variée. La grande liberté que vous m’avez laissée dans mon travail a pu être déstabilisante.
Elle m’a néanmoins fait prendre conscience de la vaste exploration qu’est la recherche
scientifique et m’a permis de mieux me connaitre et d’apprendre énormément. Merci enfin
pour m’avoir fait confiance pour mener à bien ces travaux de recherche lorsque je suis arrivé
au LEGI, avec des idées très approximatives sur la mécanique des fluides.
J’ai eu la chance de trouver une aide inestimable dans la collaboration avec Michel
Riondet pour les aspects expérimentaux. Sa connaissance des outils de mesures du LEGI
et sa disponibilité m’ont permis de mener à bien les mesures sur les modèles réduits. J’espère
pouvoir conserver et exploiter le soin et la rigueur que j’ai appris et travaillant à tes côtés.
L’ordre dans lequel apparaissent mes remerciements ne saurait être le résultat d’une
quelconque hiérarchisation de leur importance. Cette remarque me permet de remercier ici
mes parents pour tout ce qu’ils m’apportent depuis ’quelques’ années. Merci pour votre
patience et votre dévouement. Même si l’ensemble des détails de cette thèse peuvent vous
laisser perplexes, vous êtes sans aucun doute les premiers à remercier pour ma curiosité.
Vous m’avez transmis le goût pour le raisonnement et l’esprit critique dont découle ma
passion pour la recherche. Et erci Mimi pour avoir supporté toute la place que je prends
depuis vingt ans et malgré tout faire n’importe quoi quand j’en ai besoin.
Enfin, les dernières lignes seront pour tous les copains qui m’ont fait grandir au cours de
ces trois années et dont l’énumération ne sera malheureusement pas exhaustive, un grand
merci ..
Aux copains du labo pour les débats aussi interminables que récurrents, les pauses
cafés salutaires, les repas dépannés dans les moments d’oubli, les volleys, skis ou autres
barbecues. En particulier à Antoine pour m’avoir supporté dans le bureau et répondu aux
millions de questions pour lesquelles j’ai fait appel à ses lumières.
Aux nombreux habitants de la maison de Gières, de plus ou moins courte durée, pour
tout ce temps passé ensemble entre la musique, le pain, le fromage et les débats scientifiques
tableau à l’appui, les cours sur le chant des baleines, le partage des fourchettes, le réveil de
la trompette ou les exposés en terre.
Aux plus nombreux encore habitants de St-Martin-le-Vinoux, qui m’ont recueilli après
l’abandon des précédents, et m’ont offert un chouette repère aux pieds du Néron, où ont
étaient écrites quelques unes des pages de ce mémoire dans la douce musique de la meuleuse
et le braillement des brebis. Merci de nous avoir supportés, patiemment, sur la fin, et vous
être occupé de nous, voire nourris, pendant les moments critiques.
A Flo, pour ces discussions sans fin, ces chambres partagées, ta sérénité inébranlable,
et ton addiction pour la science et le piment.
A Manon, et à la belle mélodie qui me trotte dans la tête..
Table des matières
Introduction 9
5
6 TABLE DES MATIÈRES
2.3.4 Le modèle v 2 − f . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
2.4 La méthode de Couplage / maillage rotatif (Sliding mesh) . . . . . . . . . 48
2.4.1 Benchmark . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
2.4.2 Fonctionnement de la méthode Code_Saturne . . . . . . . . . . . . 50
2.4.3 Méthode de l’Université de Manchester dans Code_Saturne . . . . 53
2.5 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
2.5.1 Physique du cylindre en rotation . . . . . . . . . . . . . . . . . . . 54
2.5.2 Méthodologie de validation . . . . . . . . . . . . . . . . . . . . . . . 57
2.6 Cas test 1 : Cylindre 2D en écoulement laminaire . . . . . . . . . . . . . . 58
2.6.1 Maillage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
2.6.2 Cylindre immobile . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
2.6.3 Cylindre rotatif . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
2.7 Cas test 2 : Cylindre Turbulent . . . . . . . . . . . . . . . . . . . . . . . . 69
2.7.1 Description des simulations . . . . . . . . . . . . . . . . . . . . . . 69
2.7.2 Résultats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
2.7.3 Discussion et Comparaison avec Fluent . . . . . . . . . . . . . . . . 77
2.8 Simulation de l’hydrolienne Achard10 en 2D . . . . . . . . . . . . . . . . . 79
2.8.1 Maillages . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
2.8.2 Paramètres de calcul . . . . . . . . . . . . . . . . . . . . . . . . . . 80
2.8.3 Résultats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
2.8.4 Comparaison Code_Saturne/ Fluent / PIV . . . . . . . . . . . . . 81
Bibliographie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
Le besoin d’énergie croissant créé par une consommation toujours plus mise en avant
devient un problème important de par son impact environnemental et social. Deux grandes
façons, ‘dites durables’, d’y répondre prédominent : la conception de nouveaux systèmes de
génération électrique, sans utilisation de matières fossiles, qui sont dit renouvelables, et la
diminution de la consommation. La seconde est un processus plus lent dont les implications
sociétales ne sont pas maîtrisées, et moins avantageuses en termes de création de valeur
marchande. Cela étant dit, la combinaison de technologies propres et plus efficaces pour
l’énergie et une diminution de la consommation sont nécessaires de façon conjointe pour
répondre au défi que représente la préservation de l’environnement.
Le développement de l’énergie hydrolienne participe à cet effort en venant s’ajouter
aux dispositifs de production d’énergies renouvelables existants. Elle consiste à utiliser
l’énergie cinétique transportée par les masses d’eau en mouvement. L’énergie contenue par
un m3 d’eau qui se déplace (mesurée en joules J) est proportionnelle au carré de sa vitesse.
La puissance (en Watt W ) est proportionnelle au débit de cette énergie, c’est à dire à
la vitesse au cube. Il est donc avantageux de chercher des sites où la vitesse du courant
est la plus forte. L’implantation des hydroliennes vise donc les fleuves, mais également
les courants des marées, ou les courants océaniques, en privilégiant des sites pour lesquels
la topologie de l’environnement accélère l’écoulement (îles, rétrécissements, estuaires). Les
systèmes qui permettent la récupération de cette énergie sont typiquement de forme et
de fonctionnement semblables à ceux des éoliennes, maintenant bien connus. Bien que les
analogies soient nombreuses entre les deux technologies, on peut noter certains avantages
de l’énergie hydrolienne. L’impact visuel et le bruit produits, qui sont les principales sources
9
10 Introduction générale
des reproches faits aux éoliennes, sont quasiment absents avec les hydroliennes. De plus, la
masse volumique de l’eau plus de huit cent fois plus importante que celle de l’air fait que les
systèmes peuvent produire une énergie équivalente pour une taille trois fois inférieure. La
prédictibilité de la ressource, qui a un impact important sur la gestion du réseau électrique
est meilleure, les courants marins ou fluviaux étant plus réguliers que la vitesse du vent.
L’impact environnemental de tels systèmes est encore mal connu. Les principales études
ont jusqu’à maintenant portées sur l’effet de l’augmentation de la turbulence de l’écoule-
ment sur le lit sédimentaire, sur l’impact sur la faune et sur l’altération des courants en
cas d’extraction d’énergie trop importante.
L’effet futur sur l’environnement ne peut être nié, et une modification, au moins locale
aura lieu. Il s’agit de surveiller ses implications à moyen et long terme. On trouvera un
inventaire des études menées jusqu’à présent dans les travaux de Boehlert et Gill (2010) 1 .
L’énergie hydrolienne est souvent confondue avec les barrages marémoteurs comme celui
de la Rance, à côté de St Malo. Leurs fonctionnements sont réellement différents. Dans
le cas d’un barrage, une importante structure de génie civil est construite afin d’utiliser
l’énergie potentielle de l’élévation du niveau de l’eau dans l’estuaire. A marée montante,
l’eau rentre dans le barrage par des ouvertures, qui sont ensuite fermées au niveau le plus
haut. A marée descendante la différence de niveau crée une différence de pression qui est
ensuite turbinée. L’importante modification des courants et l’envasement de la baie ont
participé à la vision négative de ce projet. Les hydroliennes présentent l’avantage de ne
ralentir que localement le courant.
Deux grands types de turbines hydro-cinétiques existent, qui sont inspirées des éo-
liennes. Celles à flux axial et transverse. La jeunesse de la technologie et le faible retour
d’expérience fait que de nombreux designs sont actuellement au stade de prototype ou
d’essai sur les sites pilotes. Les turbines à flux transverse présentent un certain nombre de
contraintes dans l’air, surtout structurelles, qui ont limité leur utilisation. Pour les applica-
tions fluviales, ces machines avec un axe de rotation orthogonal à l’écoulement connaissent
1. Boehlert, G. W., Gill, A. B. (2010). Environmental and ecological effects of ocean renewable energy
development : a current synthesis. Oceanography, Vol.23, n2
Introduction générale 11
un nouvel essor. Le projet de recherche HARVEST, initié au début des années 2000, a
permis le développement d’une machine à axe vertical hydrodynamiquement, mécanique-
ment et électriquement fonctionnelle. A partir du concept de turbine Darrieus datant des
années 30, un nouveau design a vu le jour, sous le nom de turbine Achard. L’intégration
complète de ces recherches a débouché sur la construction d’un prototype à deux colonnes
carénées. L’application fluviale, industrialisée par la société Hydroquest, est montée sous
une barge flottante. La partie électrique peut ainsi être placée hors de l’eau, ce qui facilite
la maintenance et lui permet d’opérer dans des conditions plus favorables. Deux projets
pilotes ont permis la mise en eau de ces turbines, en Guyane et à Orléans.
La prochaine étape de ce développement est l’implantation de ces machines en ferme
fluviales ou maritimes. L’utilisation des hydroliennes requiert une vitesse de courant mi-
nimum de l’ordre de 2 m/s, avec pour autant un taux d’objet en suspension faible et une
hauteur d’eau suffisante. Les sites favorables sont peu nombreux et on essaiera d’optimiser
leur utilisation en multipliant les machines. Il est donc fondamental de comprendre l’in-
teraction des turbines et son impact sur le rendement. Le sujet des parcs a été étudié de
façon exhaustive dans le cas des éoliennes. Néanmoins, l’écoulement de couche limite at-
mosphérique et celui des rivières est différent en termes de structure (turbulence, direction,
variabilité), mais également de géométrie (confinement, hauteur d’eau, rugosité). Une ana-
lyse qualitative des effets de parc permet de comprendre l’importance de sa compréhension.
Un générateur hydro-cinétique idéal dans un écoulement en milieu infini (ie sans limites)
peut, d’après la théorie de Betz, récupérer au maximum 16/27e de la puissance du courant.
La proximité des turbines dans la direction transverse à l’écoulement aura pour effet de
forcer le débit dans les rotors et d’augmenter le rendement global. On peut ainsi définir la
première génération de parc, pour laquelle les turbines sont placées sur une ligne ou deux
lignes de manière intercalée. La seconde génération intervient lorsqu’une troisième rangée
vient placer des machines dans le sillage des premières. Elles voient alors un déficit de
vitesse et une augmentation du niveau de turbulence, dont l’effet est probablement négatif
sur le rendement et sur la tenue mécanique des structures.
Il est donc essentiel de comprendre les phénomènes en jeu dans les parcs et de disposer
12 Introduction générale
pour cela d’un outil de modélisation des turbines et de l’écoulement pour un site réel.
La complexité des géométries mise en jeu et l’écart entre les différentes échelles rendent
les représentations complètes très coûteuses, voir inenvisageables. Pour cette raison, un
modèle simplifié de la turbine est nécessaire. Il doit pouvoir reproduire les caractéristiques
complexes de la machine à partir de paramètres simples.
L’étude présentée dans ce mémoire de thèse vise la construction d’un modèle de ce
type. Elle souligne notamment la méthodologie de construction, qui passe par une étude
de la turbine réelle à l’aide d’outils de calculs complets et de mesures expérimentales. On
présentera ici les propriétés des turbines à axe vertical en général, mais également le travail
effectué sur le logiciel Code_Saturne afin de pouvoir simuler ces hydroliennes. La mise en
place du modèle pour un calcul précis de la puissance et du sillage est ensuite décrite, et son
fonctionnement validé, pour enfin estimer les effets dans des situations de parcs théoriques.
Présentons rapidement ici les différentes parties détaillées dans ce mémoire.
Le fonctionnement des turbines à flux transverse est régi par un écoulement complexe
qui a fait l’objet de nombreuses études. Le premier chapitre présente une revue des dif-
férentes évolutions du modèle original proposé par George Darrieus puis une analyse des
phénomènes physiques agissant sur l’entraînement des pales. Elle soutient l’explication des
résultats discutés par la suite. Nous présentons également les outils expérimentaux qui
permettent la caractérisation de la turbine Achard au LEGI (Laboratoire des Ecoulements
Géophysiques et Industriels).
L’étude numérique de la turbine est un outil indispensable pour la prédiction et la com-
préhension des propriétés de l’hydrolienne. Des études précédemment menées ont montré la
capacité des modèles de turbulence RANS à reproduire la dynamique tourbillonnaire dans
le rotor, avec notamment le logiciel Fluent, code général de résolution des équations de
Navier-Stokes par la méthode des volumes finis. Ces équations sont également résolues par
Code_Saturne, qui présente l’avantage d’être un outil de distribution open-source et qui est
développé par EDF. Il intègre une méthode de maillage rotatif permettant la simulation des
géométries en mouvement. Le chapitre 2 présente cette méthode et son fonctionnement. Sa
validité est évaluée par l’étude d’un cylindre en rotation dans des écoulements laminaires
Introduction générale 13
sagé dans un canal peu large. On peut ainsi déterminer l’influence des rangées successives
sur le rendement global. Cette dernière partie est également l’occasion d’évoquer les mo-
difications envisageables pour améliorer l’applicabilité du modèle.
Chapitre 1
15
16 CHAPITRE 1. LES HYDROLIENNES À FLUX TRANSVERSE
rebondissant sur le choc pétrolier, le prototype conçu et testé par les laboratoires SANDIA
(USA), a marqué le début de l’utilisation des éoliennes à axe vertical. La qualité des mesures
sur cette éolienne en soufflerie en a fait une expérience de référence pour les développements
numériques durant de nombreuses années. Par la suite, les turbines à flux transverse ont été
limitées aux machines aériennes de petite taille en raison des difficultés d’ordre mécanique
[2], et reprennent leur essor à la fin du 20e siècle avec les applications maritimes ou fluviale.
Les turbines à flux transverse (en anglais cross flow water turbine) présentent l’avan-
tage de la multidirectionnalité. C’est à dire que leur fonctionnement est indépendant de
la direction de l’écoulement, mais également du gradient de vitesse vertical présent dans
les écoulements de couche limite atmosphérique. De plus, elle présente une relative simpli-
cité pour la fabrication, notamment pour les pales, à section constantes, mais aussi pour
la maintenance, en autorisant le placement de la génératrice au niveau du sol (pour les
éoliennes) ou au niveau de la surface libre (pour des hydroliennes flottantes). Elles ont
néanmoins un rendement inférieur à celui des turbines à flux axial, mais surtout une bien
plus forte sollicitation mécanique cyclique sur la structure, et un coefficient de trainée
environ 1,5 fois plus important.
pour reprendre les forces centripètes. Les efforts sont alors transmis en tension, situation
qui offre la meilleure résistance des matériaux. Elle est testée par les laboratoires SAN-
DIA, qui obtiennent un rendement optimum de 33% environ, pour un paramètre d’avance
compris entre 4 et 5 (rapport entre la vitesse des pales et la vitesse de l’écoulement).
Figure 1.2 – Vue d’artiste de différents concepts développés depuis la Darrieus, tiré de
Islam [2].
de pale et aux raccords bras-pale, ainsi que dans la plus faible tenue mécanique pour les
designs avec un seul bras.
De l’éolien à l’hydrolien L’utilisation des turbines dans l’eau plutôt que dans l’air a
un effet sur le design et le régime de fonctionnement. La différence vient du changement de
masse volumique et de viscosité dynamique, ainsi que les plages de vitesse typiques dans
lesquelles fonctionnent ces systèmes. Le tableau 1.1, résume les paramètres pour les deux
milieux (valeurs à 20◦ C) :
Table 1.1
On trouve un rapport qualitatif de 1/50 pour les forces centrifuges pour des machines
de même taille, ce qui explique que ces forces ne sont plus le paramètre dominant de
dimensionnement (comme il l’est pour une forme troposkienne) qui considérant plutôt
l’effet du chargement en pression et notamment les chargements instationnaires.
Les évolutions majeures Les évolutions majeures apportées au concept Darrieus d’ori-
gine portent essentiellement sur la périodicité du couple. Une pale au cours de la rotation
voit alternativement des zones de couple moteur et frein. Lors de l’utilisation nominale, la
somme des trois contributions ainsi que l’inertie de la machine à tendance à lisser le couple
récupéré.
La turbine Gorlov [4] a été conçue en prenant en compte les inconvénients de la Darrieus
cité précédemment. Afin de répartir les efforts, il impose une forme hélicoïdale à des pales
verticales (c’est à dire de rayon constant), comme illustré sur la figure 1.3. La projection
de la turbine sur un plan horizontal montre ainsi un recouvrement plus grand du cercle
de rotation. Le recouvrement est complet dans le cas d’une machine seule, et partiel pour
une tour de turbine, pour lequel les vides sont comblés par un déphasage des différentes
machines sur la tour. Ce concept est attirant mais les résultats en termes de performance
sont assez approximatifs. L’étude de Shiono [5] montre une dégradation du coefficient de
puissance CP jusqu’à 18%, par rapport aux 35% avancés par Gorlov.
La turbine à pale droite permet de modifier l’angle de calage φ des pales (angle entre
la directrice du profil et la tangente au cercle de rotation). De nombreux systèmes ont été
proposés pour ajuster l’incidence des pales durant la rotation (ils sont dits à pitch variable).
Le mouvement peut être imposé par le système en fonction de la position angulaire, il
nécessite dans ce cas un apport d’énergie et un système de contrôle sophistiqué. L’option
peut également être prise de laisser un degré de liberté en rotation autour d’un axe, pour une
plage de φ limitée, typiquement φ = ±3◦ . Le moment sur le profil contribue à augmenter
sa plage motrice. Le mouvement est destiné à limiter l’angle d’incidence du profil et par ce
biais l’apparition du décollement [7]. Le rendement peut ainsi être augmenté bien que cette
amélioration théorique n’est pas toujours constatée. On remarque néanmoins une capacité
d’auto-démarrage qui n’est pas forcément présente pour les pitchs fixes. La complexité
supplémentaire induite par de tels systèmes met en péril l’objectif originel de simplicité de
20 CHAPITRE 1. LES HYDROLIENNES À FLUX TRANSVERSE
Figure 1.3 – (a) Turbine Gorlov, tiré de Gorban [6], (b) tour bi-rotor de Hydroquest
la turbine Darrieus.
façon générale, le principe physique reste inchangé. Cette partie décrit le fonctionnement
général des turbines cross-flow. La description est faite en 2D dans un plan orthogonal à
l’axe de rotation.
Fl
α
Bord
d'attaque Extra Bord de
Epaisseur dos fuite
Ligne
Cambrure Moyenne
Corde
Intrados
a) b)
Figure 1.4 – a) Terminologie d’un profil aérodynamique, (b) Schéma des forces de trainée
et portance.
Le principe des turbines à flux transverse consiste à récupérer l’énergie induite par les
forces qui s’exercent sur un profil en rotation. D’un point de vue général, un profil de
pale a une forme aérodynamique optimisée pour les efforts qu’elle subit dans un fluide en
mouvement. Ces efforts sont dus à la déviation que le profil impose à l’écoulement. La
figure 1.4.a) reprend la terminologie utilisée en aérodynamique. Le profil est une forme
longiligne avec le point amont dans le fluide, appelé bord d’attaque, et le point aval, le bord
de fuite. Ces deux points sont reliés par la ligne moyenne. La cambrure du profil est l’écart
entre la ligne moyenne et la ligne droite reliant les deux extrémités (appelé corde). Elle est
caractéristique de l’asymétrie du profil. La face interne du profil est appelée intrados, et la
face externe extrados. L’épaisseur du profil est définie comme la largeur du profil dans la
direction orthogonale à la ligne moyenne. La position du maximum de l’épaisseur joue un
rôle important dans la transition de couche limite laminaire à turbulente.
Pour un profil donné, la force exercée dépend de l’angle d’incidence et du régime d’écou-
lement. L’angle d’incidence, noté α est défini comme l’angle entre la corde du profil et la
−
→
direction du fluide incident W dans le référentiel de la pale. α est positif quand l’écoulement
voit l’intrados.
→
−
La force est usuellement projetée en deux composantes : la force de portance Fl (lift),
−
→
orthogonale à la vitesse W , et prise positive en s’éloignant de l’extrados, et la force de
22 CHAPITRE 1. LES HYDROLIENNES À FLUX TRANSVERSE
1.4
1.2
0.8
0.6
0.4
0.2
0
0 5 10 15 20 25 30
a) b)
Figure 1.5 – (a) Régimes d’écoulement sur un profil en écoulement turbulent, (b) In-
fluence du nombre de Reynolds sur les coefficients de portance (lignes pleines) et de trainée
(pointillées), d’après les mesure de Sheldahl [9]
−
→ −
→
trainée Fd (drag), dans la direction et le sens de W . Dans le cas d’un aéronef, muni d’ailes
profilées, la portance s’oppose à la chute par pesanteur, et la trainée ralentit son mouve-
ment. Le rapport Fl /Fd est appelé finesse du profil. La figure 1.5.b) présente les coefficients
de portance et de trainée pour un profil NACA0018 à différents nombre de Reynolds. On
appelle coefficient aérodynamique l’adimensionnalisation des forces par une force de ré-
férence, en général la pression dynamique qui s’applique sur une surface équivalente à la
corde c multipliée par l’envergure de la pale L, on définit CL et CD par :
FD,L
CD,L = 1 (1.3)
2
ρcLU02
La portance augmente avec l’incidence α jusqu’à une valeur critique, appelé angle de décro-
chage. Jusqu’à cette valeur, la couche limite est attachée à l’extrados. Au delà, on constate
l’apparition d’un point de décollement causé par un gradient de pression adverse à la sur-
face. Le coefficient de trainée montre d’abord un palier jusqu’à la valeur de décrochage,
à partir de laquelle il augmente rapidement, suivant le déplacement du point de décolle-
ment du bord de fuite vers le bord d’attaque. Une fois la couche limite totalement décollée,
l’augmentation continue, bien que plus douce. Le coefficient de portance augmente jusqu’au
décrochage puis chute.
1.2. ANALYSE PHYSIQUE DES TURBINES À FLUX TRANSVERSE 23
-V W
D U0
L
Ft,res
U0
D
-V R
Ft,motrice
W
L 1
Y
2 X
D
Ft,res
L
W
Ft,motrice
-V
U0
Figure 1.6 – Vue schématique de la direction des efforts sur une pale en rotation
néfastes à la puissance produite. Au-delà de l’optimum, la portance est limitée et les ef-
fets visqueux deviennent prépondérants. Les frottements limitent alors le rendement de la
machine. On distingue donc 3 zones que l’on retrouve sur la Figure 1.7 b) : la zone de
décrochage dynamique, l’optimum, et la zone des effets secondaires.
Ces trois comportements de la dynamique de l’écoulement présentent chacun des défis
différents pour la simulation. La zone de décrochage dynamique présente des écoulements
fortement décollés que les modèles RANS ont des difficultés à reproduire avec précision.
On obtient de meilleurs résultats pour les simulations réalisées autour de λ =2-3, valeurs
pour lesquelles les vortex sont moins développés.
35
nsition
50
30 Tra
Eff
40
et
iqu
30 25
ss
m
ec
20
na
CP (%)
20
on
dy
10
da
ge
ire
0 15
ha
s
oc
−10
cr
10
Dé
−20
−30 5
−40
0
−50
0 0.5 1 1.5 2 2.5 3
λ
0 50 100 150 200 250 300 350
a) b)
décrochage est plus grand, et la portance est augmentée. La chute de portance et l’aug-
mentation de trainée associée est plus brutale. La variabilité du chargement et donc de la
fatigue sur les pales est d’autant augmentée.
Nous donnerons ici une description du décrochage dynamique et l’influence qu’il a sur
une turbine à flux axial. Le lecteur pourra trouver une revue approfondie de la littérature
dans la thèse de Jonathan Bossard [10] qui a étudié le décrochage dynamique sur une
turbine Darrieus par la technique de vélocimétrie par laser doppler (PIV).
πf c
kr = (1.9)
2V0
Ce second paramètre permet d’opérer une distinction entre les deux régimes du décrochage
dynamique. Le décrochage dynamique est dit léger lorsque αmax est à peu près égal à
αstatique , et dit profond, quand il le dépasse largement.
Le décrochage dynamique léger est similaire au décrochage classique. Comparativement
au cas statique, la portance est augmentée dans la partie α croissant qui précède le décro-
chage. La partie décroissante de α est caractérisée par un CL plus faible qu’à la montée.
Sur la figure 1.8 pour un profil NACA0012, cet hystérésis est marqué entre 3◦ et 15◦ . Mc-
Croskey et al. [11] identifient deux comportements de l’écoulement pour ce régime. Pour
un profil moyen, le décollement intervient au bord d’attaque et se propage sur toute la
surface, pour les profils plus épais, le décollement apparaît classiquement au bord de fuite,
et le point de décollement migre vers le bord d’attaque. Le premier cas est caractérisé par
une variation plus brutale des efforts.
1.8 0.4
αup
αdown
αstat 0.3
0.8
0.2
CD
CL
0.1
−0.2
0
−1.2 −0.1
−15 −10 −5 0 5 10 15 20 −15 −10 −5 0 5 10 15 20
α (°) α (°)
Figure 1.8 – Illustration du décrochage léger sur un profil oscillant par les coefficients de
portance et trainée. D’après Lee et Gerontakos [12]
Application aux turbines Darrieus Dans les turbines Darrieus, les conditions sont
requises pour observer un décrochage dynamique. Le phénomène est globalement compa-
rable au profil oscillant. Néanmoins, les effets tridimensionnels peuvent être très impor-
tants. Les parties marginales (extrémités de pale, bras) engendrent des vitesses verticales
qui déclenchent le décollement. D’après Zanette [14], la fréquence critique pour laquelle le
mouvement ne peut plus être considérée quasi-statique est dépassé pour toutes les solidités
supérieures à 0,05 1 . Dans le cas des applications hydrauliques, elle sera toujours dépas-
sée. Le décrochage dynamique dépend donc surtout de l’angle d’incidence maximal atteint
pendant la rotation qui est directement relié à λ. Laneville et Vittecoq [15] comparent les
coefficients de portance et de trainée pour une éolienne de solidité σ = C/R = 0, 2 à diffé-
rents λ. Le décrochage dynamique est constaté en dessous de la valeur limite λ =3. L’effet
d’hystérésis est en forte augmentation lorsque λ diminue. Les visualisations de Brochier
et al. [16] leur permettent d’illustrer la dynamique tourbillonnaire dans le sillage dans le
rotor à λ =2,14 (figure 1.10). On remarque la formation d’un tourbillon à l’extrados (in-
térieur du cercle) aux alentours de θ = 70◦ , qui se détache vers θ = 120 − 130◦ . Ferreira
et al. effectuent des mesures similaires avec un système PIV. Cette instrumentation leur
1. La solidité de la turbine représente sa porosité. On peut montrer que la vitesse de rotation diminue
avec la solidité.
28 CHAPITRE 1. LES HYDROLIENNES À FLUX TRANSVERSE
2.5 αup ➃
αdown
2 αstat
➁➂
1.5 ➅ (a)
➀ ➄ TBL
CL
1
0.5 ➆ Rev
ers
e
ow
0 (b)
➇
−0.5
−10 −5 0 5 10 15 20 25 30
1
➃ (c)
LEV
0.8
0.6
CD
0.4 ➁
(d)
LEV
0.2 ➀
−10 −5 0 5 10 15 20 25 30
(e)
0.01
➁
0
➇ ➀➆
−0.01 (f)
CM
−0.02
−0.03
−0.04 ➃
(g)
−10 −5 0 5 10 15 20 25 30
α (°)
permet de calculer l’intensité de la circulation autour des pales. Ils montrent que la force
du tourbillon augmente tant qu’il est attaché à la pale et diminue après le décrochage. La
thèse de Jonathan Bossard a montré que les effets 3D jouent un rôle de déclencheur du
décollement en comparant les mesures PIV à des calculs 2D et 3D. Il constate l’absence de
tourbillon de bord d’attaque à partir de λ = 2, 5, ainsi qu’un second lâcher tourbillonnaire
vers θ = 120 − 130◦ .
La veine d’essai
La veine d’essai permettant l’étude des maquettes d’hydrolienne a été conçue lors du
DRT de Linda Guitet [17]. Le système de contrôle et de mesure de la turbine a été mis au
point et calibré par Nicolas Dellinger [18] et Vivien Aumelas [19].
La veine d’essai qui suit le convergent, est d’une longueur de 1000mm, et a une section
rectangulaire de largeur 700mm et de hauteur 250mm. Les faces latérales et inférieures
1.3. PROGRAMMES D’ESSAIS EXPÉRIMENTAUX 31
sont des vitres en plexiglass transparent d’épaisseur 50mm, montées sur une structure
métallique. La face supérieure est une plaque en aluminium fixée à l’armature principale.
Cette section de 0, 175m2 permet une plage de vitesse allant de 0,8 à 2,9 m/s. L’aluminium
choisi pour pouvoir déplacer la plaque facilement (à l’aide de palans fixés à la structure du
bâtiment) en raison de sa légèreté pose néanmoins un problème de rigidité. On constate
une déformation lors de la mise en eau du circuit, ce qui est problématique pour la mesure
des efforts [19].
L’instrumentation
La balance de mesures, positionnée sur la plaque supérieure assure deux fonctions :
l’asservissement de la turbine en vitesse et la mesure des efforts appliqués par le fluide.
Les éléments constitutifs sont visibles sur la figure 1.12.b). Une plaque métallique carrée
de 300mm x 300mm repose sur quatre capteurs piezo-électriques en contact avec la plaque
supérieure de la veine. L’arbre de la turbine représente l’axe de symétrie de tout le système.
Il traverse les différents étages et notamment la plaque supérieure. La turbine est fixée à
l’arbre par le moyeu à mi-hauteur de la veine. Il est également fixé au rotor de la génératrice
et du résolveur. Le guidage en position est assuré par le boîtier à roulement, dont le bâti
repose sur la plaque carrée. La génératrice permet d’imposer une vitesse de rotation à
la turbine (fonctionnement en moteur) et de mesurer le couple récupéré. Le corps de la
32 CHAPITRE 1. LES HYDROLIENNES À FLUX TRANSVERSE
∂Ω
J = Ct + Cg + Cf (1.10)
∂t
Le couple de la génératrice est directement proportionnel à l’intensité mesurée. Le
couple de frottement est proportionnel à la vitesse de rotation mais présente une faible
répétabilité. Il est évalué par mesure du couple en l’absence de turbine avec une précision
de l’ordre de 2%. L’inertie de la machine participe à lisser le couple reçu par la turbine.
Elle est en moyenne nulle sur un tour en régime établi, mais doit être prise en compte pour
évaluer le couple instantané. L’incertitude globale sur la mesure du couple est estimée à
4% [19].
Figure 1.12
a) A 10 b) A 12
40
35
30
25
20
15
10
0
0 0.5 1 1.5 2 2.5 3
Figure 1.13 – Design et performance des turbines A10 et A12 (d’après Bossard [10])
1.3. PROGRAMMES D’ESSAIS EXPÉRIMENTAUX 35
namique de l’écoulement. Les maquettes sont étudiées dans une plage de fonctionnement
comprise entre λ = 1 et λ = 4. A l’optimum, le paramètre d’avance de 2 équivaut à une
fréquence de rotation de 8,6Hz pour U∞ = 2, 3 m/s.
λU∞ c
ReC = ρ (1.11)
µ
B
A
produit également une amélioration du rendement [23]. Dans ce cas le plan de symétrie
entre deux turbines joue le rôle d’une paroi qui bloque le fluide. L’effet de contournement
rentre alors en jeu à l’échelle de la ligne de turbine.
Pour étudier cet effet, un dispositif expérimental a été conçu au cours de cette thèse. Un
jeu de plaques verticales (F) de longueur 625 mm est placé dans le sens de l’écoulement de
manière symétrique de part et d’autre du modèle réduit (voir figure 1.14). Elles simulent un
rétrécissement de la veine. L’extrémité amont des plaques se trouvent 130 mm avant l’axe
de rotation. Le bord d’attaque présente un rayon de courbure de 7,5 mm. Les simulations
numériques préliminaires réalisées ont montré l’absence de recirculation sur ces parois. La
turbine crée une perte de charge dans la partie centrale. Afin de forcer le flux, un jeu de
plaques à trous (E) est placé orthogonalement au fluide dans la partie du contournement.
Ce système est fixé dans une plaque de support supérieure horizontale (A et B) conçu par
Ane Mentxaca pour l’adaptation des carénages. Il est plaqué sur la vitre inférieure, avec
une bande de mousse antidérapante.
Le débit dans la section centrale est mesuré par laser LDV. La vitesse moyenne dans le
rotor est également mesurée, ce qui permet de connaître la vitesse de rotation spécifique
globale et locale (ie adimensionnalisée par la vitesse interne). Le couplage de ces mesures
avec celles des efforts ne peut être présenté ici par faute de temps mais fait actuellement
l’objet d’un travail visant à valider les résultats des simulations RANS de la turbine confi-
née.
BIBLIOGRAPHIE 37
Bibliographie
[1] G. Darrieus, “Turbine having its rotating shaft transverse to the flow of the current,”
Dec. 8 1931. US Patent 1,835,018.
[2] M. Islam, D. S.-K. Ting, and A. Fartaj, “Aerodynamic models for darrieus-type
straight-bladed vertical axis wind turbines,” Renewable and Sustainable Energy Re-
views, vol. 12, no. 4, pp. 1087–1109, 2008.
[4] A. M. Gorlov, “Unidirectional helical reaction turbine operable under reversible fluid
flow for power systems,” Sept. 19 1995. US Patent 5,451,137.
[5] M. Shiono, K. Suzuki, S. Kiho, et al., “Output characteristics of darrieus water turbine
with helical blades for tidal current generations,” in Proceedings of the Twelfth Inter-
national Offshore and Polar Engineering Conference, Kitakyushu, Japan, pp. 859–864,
2002.
[6] A. N. Gorban, A. M. Gorlov, and V. M. Silantyev, “Limits of the turbine efficiency for
free fluid flow,” Journal of Energy Resources Technology, vol. 123, no. 4, pp. 311–317,
2001.
[7] B. Kirke and L. Lazauskas, “Limitations of fixed pitch darrieus hydrokinetic turbines
and the challenge of variable pitch,” Renewable Energy, vol. 36, no. 3, pp. 893–897,
2011.
[11] W. McCroskey, “The phenomenon of dynamic stall.,” tech. rep., DTIC Document,
1981.
[12] T. Lee and P. Gerontakos, “Investigation of flow over an oscillating airfoil,” Journal of
Fluid Mechanics, vol. 512, pp. 313–341, 2004.
[15] A. Laneville and P. Vittecoq, “Dynamic stall : the case of the vertical axis wind
turbine,” Journal of Solar Energy Engineering, vol. 108, no. 2, pp. 140–145, 1986.
[19] V. Aumelas, Modélisation des hydroliennes à axe vertical libres ou carénées : développe-
ment d’un moyen expérimental et d’un moyen numérique pour l’étude de la cavitation.
PhD thesis, Grenoble-INP, 2011.
[21] A. Bahaj, A. Molland, J. Chaplin, and W. Batten, “Power and thrust measurements
of marine current turbines under various hydrodynamic flow conditions in a cavitation
tunnel and a towing tank,” Renewable energy, vol. 32, no. 3, pp. 407–426, 2007.
[22] A. Mentxaca, “Analyse numérique des hydroliennes à axe vertical munies d’un caré-
nage,” 2011.
[23] S. Antheaume, T. Maître, and J.-L. Achard, “Hydraulic darrieus turbines efficiency for
free fluid flow conditions versus power farms conditions,” Renewable Energy, vol. 33,
no. 10, pp. 2186–2198, 2008.
Chapitre 2
39
40 CHAPITRE 2. LE MAILLAGE ROTATIF DANS CODE_SATURNE
Code_Saturne développé par EDF. Une collaboration étroite avec les équipes de dévelop-
pement du département R&D de EDF, a permis de mettre en évidence les capacités du
solveur Code_Saturne.
L’objectif de ce chapitre est de présenter Code_Saturne, et la méthode de maillage
glissant qui permet de simuler l’écoulement autour de l’hydrolienne en rotation. Le cas du
cylindre en rotation dans des écoulements laminaire et turbulent, abondamment étudiés et
documentés, a été utilisé pour la validation de la méthode. Celle-ci a été ensuite appliquée
à l’étude de l’hydrolienne Achard en approximation bidimensionnelle dont la configuration
correspond à celle de la veine d’essai expérimentale du LEGI. L’objectif était de mener
une étude extensive du comportement des hydroliennes en 2-D. Il était nécessaire d’étudier
l’influence du comportement du confinement dans le tunnel ainsi que l’impact du sillage
sur l’écoulement.
procédure. Pour une variable donnée, un critère d’arrêt est fixé sur la différence relative
entre deux itérations successives. Les étapes de l’algorithme de résolution sont détaillées
ci-dessous :
• L’étendue des échelles : la turbulence est présente sur un large spectre d’échelles
spatio-temporelles, depuis des tourbillons de la taille du domaine disponible pour
l’écoulement, jusqu’à la taille minimale à laquelle les tourbillons disparaissent en
agitation moléculaire, traduite au niveau macroscopique par la viscosité du fluide.
La distance entre les extrémités du spectre augmente avec le nombre de Reynolds de
l’écoulement.
• La cohérence des structures. Les structures turbulentes (ou tourbillons) sont obser-
vables et les corrélations entre les composantes des fluctuations de vitesse sont non
nulles.
L’énergie du champ moyen de vitesse est dissipée par la turbulence. Les plus grands tour-
billons (qui transportent le plus d’énergie) se brisent et forment des tourbillons plus petits
et ainsi de suite (cascade des échelles), jusqu’à ce que les effets visqueux soient prépondé-
rants. A cette échelle, l’énergie est dissipée sous forme de chaleur. Cette cascade d’énergie
des grandes échelles vers les petites, porte le nom de Kolmogorov qui l’a expliqué en 1942
[12].
u = U + u′ avec u = U et u′ = 0 (2.3)
D’un point de vue physique, la moitié de la somme des variances des composantes de la
vitesse fluctuante est interprétée comme l’énergie cinétique contenue dans la turbulence, ou
2.3. LES MODÈLES RANS 43
1
k = (u′2 + v ′2 + w′2 ) (2.4)
2
et
q
1
( 23 k)1/2 3
(u′2 + v ′2 + w′2 )
It = = (2.5)
Uref Uref
On réécrit les équations de conservation de l’énergie avec cette décomposition puis on
applique l’opérateur de moyenne pour faire disparaître les termes fluctuants de premier
ordre. Par (2.1) on obtient :
∂Ui ∂Ui ′
′ ∂ui ∂P ∂ 2 Ui
ρ[ + Uj + uj ] = ρFi − +µ (2.7)
∂t ∂xj ∂xj ∂xi ∂xj ∂xj
∂u′
Par application de (2.6), on a u′j ∂xji = ∂x∂ i u′j u′i . Le tenseur défini par u′j u′i est appelé tenseur
de Reynolds. Les termes diagonaux de ce tenseur sont des contraintes normales non-nulles,
car égaux au carré des fluctuations de vitesse. Les autres termes sont des contraintes de
cisaillement, ils découlent de la corrélation entre les composantes de la vitesse. Ils sont éga-
lement non nuls par la nature cohérente des tourbillons. On notera que ces contraintes sont
généralement supérieures en intensité aux contraintes visqueuses en écoulement turbulent.
La non-linéarité des termes d’advection dans l’équation de départ implique l’obtention
d’un système d’équations ouvert (plus d’inconnues que d’équations) pour lesquelles on va
devoir trouver des relations de fermeture. On peut ici distinguer différentes solutions pour
la fermeture de ce système :
• Les modèles du premier ordre, dans lesquels l’effet des tensions de Reynolds est
modélisé à l’aide d’une longueur de mélange, fixe ou dépendante de l’écoulement.
• Les modèles du second ordre, où les tensions de Reynolds sont calculées directement,
les moments d’ordre supérieur sont modélisés.
44 CHAPITRE 2. LE MAILLAGE ROTATIF DANS CODE_SATURNE
Nous poursuivrons l’analyse en nous concentrant sur les modèles du premier ordre.
Boussinesq propose en 1877 l’hypothèse que les contraintes turbulentes sont proportion-
nelles aux vitesses de déformation du champ moyen. La proportionnalité fait intervenir µt ,
la viscosité turbulente.
∂Uj ∂Ui 2
−ρu′j u′i = µt ( + ) − ρkδij (2.8)
∂xi ∂xj 3
Cette hypothèse est erronée en ce qu’elle suppose une isotropie de la turbulence. On définit
S le tenseur des taux de déformations basé sur le champ moyen RANS, qui est défini par
ses composantes :
1 ∂Uj ∂Ui
Sij = ( + ) (2.9)
2 ∂xi ∂xj
et s le tenseur des déformations basé sur la partie fluctuante de la vitesse :
∂u′j ∂u′i
s′ij =( + ) (2.10)
∂xi ∂xj
Cette hypothèse permet de linéariser les équations sur la vitesse. Le problème de ferme-
ture consiste maintenant à calculer une solution pour la viscosité turbulente. On cherche
une équation régissant l’évolution de l’énergie cinétique turbulente. En prenant le produit
scalaire de la vitesse fluctuante u′ avec l’équation de Navier-Stokes, on obtient l’équation
suivante :
∂k ∂k ∂ 1
ρ[ + Ui ]= (− p′ u′j + 2µu′i s′ij − ρu′i u′i u′j ) − 2µs′ij · s′ij − ρu′i u′j · Sij (2.11)
∂t ∂xi ∂xj |{z} | {z } |2 {z } | {z } | {z }
P ression visqueuses Reynolds Dissipation P roduction
2.3.1 Le modèle k − ǫ
Le modèle k − ǫ est le plus ancien et le plus répandu de ces modèles [13]. Il reprend
l’équation exacte développée pour k dans laquelle les termes de transport sont modélisés
avec une hypothèse de transport par diffusion de gradient. Il propose ensuite une équation
de conservation pour ǫ construite en miroir par rapport à celle sur k. Les deux équations
de conservation du modèle sont :
∂ρk µt
+ ∇ · (ρkU ) = ∇ · [ ∇k] + P − ρǫ (2.13)
∂t σk
∂ρǫ µt ǫ ǫ2
+ ∇ · (ρǫU ) = ∇ · [ ∇ǫ] + Cǫ1 P − Cǫ2 ρ (2.14)
∂t σǫ k k
avec :
P = 2µt Sij · Sij (2.15)
Le terme de production vient directement de l’application de l’hypothèse de Boussinesq au
terme P de l’équation exacte. La viscosité turbulente qui sert à fermer l’équation sur u est
déterminée à l’aide d’une analyse dimensionnelle :
k2
µt = ρCµ (2.16)
ǫ
Les constantes du modèle telles que données par Jones et Launder, et utilisées dans
Code_Saturne sont :
∂ρk µt
+ ∇ · (ρkU ) = ∇ · [ ∇k] + P − βρkω (2.19)
∂t σk
∂ρω µt ω ω2
+ ∇ · (ρωU ) = ∇ · [ ∇ω] + γ1 P − β1 ρ (2.20)
∂t σω k k
Cette modification est appelée Shear Stress Transport, dans le sens où elle imite un des
effets remarqué du terme de transport de τ dans d’autres modèle. Le modèle kω − SST
classique intègre les deux modifications.
2.3.4 Le modèle v 2 − f
Le modèle v 2 − f a été proposé par Durbin [16], afin de considérer l’anisotropie du
tenseur de Reynolds près des parois. Il permet ainsi de prendre en compte de façon plus
réaliste les écoulements détachés et impactant. La modélisation proposée rajoute deux
variables au modèle k − ǫ qui sont v 2 , la vitesse transverse à la paroi, et son terme source
f , en plus de k, l’énergie turbulente et ǫ, la dissipation. Le modèle v 2 − f est basé sur trois
équations de transport pour v 2 , k, ǫ et une équation de relaxation elliptique pour f .
Décrivons plus en détail le fonctionnement du modèle de départ proposé par Durbin en
1991. Dans les zones éloignées de la paroi, les équations de transport classiques peuvent
être utilisées. Les équations de transport pour k et ǫ sont les mêmes que dans le modèle
classique (eqs (2.13) et (2.14)). L’équation sur v 2 est définie comme suit :
Dv 2 ∂v 2 ǫ µt
= + Uj ∇j v 2 = kf − v 2 + ∇j [(µ + )∇j v 2 ] (2.26)
Dt ∂t k σk
L’équation de relaxation elliptique est résolue pour f , pour laquelle T est un temps carac-
téristique de la dissipation et L une longueur de mélange caractéristique :
1 v2 Pk
L2 ∇ 2 f − f = (C1 − 1)[ − 2/3] − C2 (2.27)
T k k
La viscosité turbulente adopte la nouvelle forme par rapport au modèle k − ǫ :
ν t = Cµ v 2 T (2.28)
2.4.1 Benchmark
La méthode Chimera est une méthode de simplification de maillage qui a commencée à
être utilisée dans les années 1980 par Benek et al [19], et Steger [20]. Elle consiste à utiliser
un maillage différent pour chacun des éléments du domaine à calculer. Considérons le cas
avec deux maillages, le premier est fixe, le second maillage est placé au dessus du premier
et peut être déplacé de façon arbitraire. Les cellules du premier maillage complètement
recouvertes sont éliminées pour le pas de temps courant, par un méthode appelée hole cut-
ting. On pourra voir les travaux de Wang [21] sur cette méthode particulière. Le processus
de construction du domaine à chaque itération est coûteux en temps de calcul, notamment
en 3D. Le couplage est assuré par la continuité des flux entre les différentes parties. Cette
méthode est très utilisée pour la simulation de géométries complexes en mouvement telle
que les rotors d’hélicoptère sur lesquels on observe des mouvements de rotation ainsi que
des battements de déformation verticale. Elle devient difficile à mettre en œuvre lorsque
l’intervalle entre rotor et stator devient petit dans le cas, par exemple, des turbomachines
ou des hélices et hydroliennes carénées.
2.4. LA MÉTHODE DE COUPLAGE / MAILLAGE ROTATIF (SLIDING MESH) 49
Maillage glissant Dans un cas plus simple, pour lequel on peut définir le mouvement
du maillage par une révolution ou une translation uniforme, la simplicité de la méthode
par maillage glissant est en général préférée. Elle utilise plusieurs maillages recouvrant des
parties distinctes de l’espace. Cette méthode est généralement utilisée pour des simulations
de profils oscillants, de turbo-machines, pompes, ou générateurs hydrocinétique. Les efforts
de développement pour de telles méthodes sont concentrés sur les points suivants :
• la parallélisation de la méthode
La principale différence entre les exemples de maillage glissant proposés réside dans
l’interpolation des variables à l’interface. Cette interpolation est faite par une méthode de
construction de cellules fantômes (“halo”) du côté opposé à l’interface. Steijl et Barakos [22]
présentent une méthode de simulation appliquée à un rotor d’hélicoptère, avec un maillage
structuré par blocs. Cela leur permet de créer deux couches de cellules ’halo’ pour le do-
maine 1 dans le domaine 2. Les variables au centre de ces nouvelles cellules sont déterminées
par interpolation dans le domaine 2. Cette technique est réservée à des maillages structu-
rés, et nécessite un stockage de données supplémentaires. Blades et Marcum [23] proposent
une méthode pour maillage non-structuré avec une couche de mailles ’halo’. Le centre de
ces mailles est placé à une distance symétrique du centre de la dernière couche de mailles
réelles par rapport à l’interface de couplage. Dans leur cas les grandeurs sont extrapolées
depuis le centre le plus proche, puis l’écoulement est résolu sur les mailles contenant les
nouvelles cellules, comme faisant partie du maillage. Elles imposent les conditions limites
pour chacun des côtés de l’interface.
Une méthode alternative présentée par Petit et. al. [24] revient à effectuer un recollement
de deux maillage à chaque pas de temps. Cette méthode dite de patching, est implémentée
dans OpenFOAM. Elle existe également dans le Code_Saturne depuis la version 3.2 (Dé-
cembre 2013). Elle assure la conservativité par construction mais est plus lourde en termes
de calcul. Comparativement, certaines méthodes de maillage glissant avec interpolation
demandent au moins deux itérations sur la résolution pour assurer la conservation de la
masse.
Immersed boundary method Peskin [25] et plus récemment Jendoubi et. al. [26] pro-
posent une méthode avec laquelle un objet fictif peut être déplacé dans un espace maillé
50 CHAPITRE 2. LE MAILLAGE ROTATIF DANS CODE_SATURNE
comme le montre la figure 2.1. A chaque itération, la géométrie de l’objet est définie en
référentiel lagrangien. Les forces appliquées sur l’objet sont calculées dans ce référentiel
puis imposées dans le champ eulérien pour les mailles qui sont intersectées par le contour.
Cette méthode initialement appliquée aux écoulements vasculaires, présente néanmoins le
problème de ne pas pouvoir construire un maillage adapté aux calculs d’écoulement décollés
en couche limite.
Figure 2.1 – Champ de vitesse dans les référentiels eulériens et lagrangiens pour la mé-
thode IBM [27]
Figure 2.2 – Schéma d’une face interne (gauche) et d’une face de bord (droite) dans le
fonctionnement classique du Code_Saturne. Extrait de [28]
• obtenir les mêmes flux convectifs et diffusifs pour une face couplée que ce qu’on aurait
pour une face interne équivalente
C2
I C3
F2
K
F3
C1
J
F1 L
F4 div(v)=0
C4
div(v)=0
div(v)≠0
div(v)≠0
La première méthode pouvant être instable, on effectue un test de pente qui permet
aI − aJ
de basculer vers la deuxième si le rapport devient trop grand. Afin de garantir la
IJ
symétrie des flux, on impose une condition limite de type Dirichlet. Cette condition est
évaluée à l’aide d’un schéma centré pour tous les scalaires, la pression et la vitesse et un
schéma upwind pour les grandeurs turbulentes.
Dans le cas général, ces conditions aux limites servent à réaliser une première étape de
prédiction de vitesse. Une étape de correction vient ensuite : les valeurs de flux de masse et
de pression au centre des cellules sont corrigées afin de respecter la divergence nulle de la
vitesse dans chaque cellule. Les valeurs de la vitesse sur chaque face sont donc modifiées.
Il en va de même pour l’interface de couplage. Or ici, la valeur de la vitesse en F est
corrigée de chaque côté de façon indépendante. Cette opération modifie la condition limite
imposée à l’étape précédente. On perd alors, a priori, l’égalité des flux. Il est donc choisi de
désactiver la correction de pression pour les mailles de bord afin d’imposer de façon exacte
la conservation du flux. Il en découle alors un calcul non conservatif.
Couplage instationnaire non-conforme Dans le cas général, les deux maillages sont
en mouvement relatif, l’un par rapport à l’autre. A chaque pas de temps, une nouvelle
correspondance entre les mailles en vis à vis est à déterminer. Comme on le voit sur le
schéma, Figure 2.3, le flux sortant par la face F 1 devrait être reporté dans les cellules 3 et 4.
Dans l’implémentation actuelle, l’interpolation (i.e. pondération du flux envoyé dans chaque
2.4. LA MÉTHODE DE COUPLAGE / MAILLAGE ROTATIF (SLIDING MESH) 53
face) n’est pas faite, la cellule 1 "verra" donc seulement la cellule 4. Cette particularité
impose donc des restrictions sur la taille des mailles dans la direction longitudinale à
l’interface. L’interpolation n’est pas faite pour ne pas alourdir le temps de calcul.
RANS et LES montrent que le déplacement du maillage par pas de temps doit rester
petit rapport à la taille des mailles des l’interface. De plus, la correction de pression est
54 CHAPITRE 2. LE MAILLAGE ROTATIF DANS CODE_SATURNE
2.5 Validation
2.5.1 Physique du cylindre en rotation
L’étude de la méthode de maillage rotatif à interface glissante impose une géométrie
d’étude à symétrie de révolution. On est amené à s’intéresser au cas du cylindre en rotation
qui remplit cette contrainte, reste géométriquement très simple et permet cependant de
conduire à des écoulements complexes. Les propriétés de ces écoulements ont été, et restent,
abondamment étudiées grâce aux moyens expérimentaux et numériques. Des compilations
de résultats faites notamment par Zdravkovitch [31] ou Roshko [32], ont synthétisé et
classifié les nombreux phénomènes en différents régimes d’écoulement.
Dans la très complète étude réalisée par Zdravkovitch, on trouve une interprétation
générique qui partitionne l’écoulement en 4 zones représentées sur la Figure 2.5 (a) :
(a)
Figure 2.5 – Zones de l’écoulement autour d’un cylindre [31] (gauche), Schémas des dif-
férents régimes d’écoulement autour d’un cylindre (droite)
2.5. VALIDATION 55
La couche limite développée est constituée de deux zones : une première, la plus éten-
due, possède un gradient positif de pression et une seconde, beaucoup plus courte, où ce
gradient devient négatif juste avant le point de détachement. Les couches limites détachées
continuent d’évoluer à l’aval comme des zones de cisaillement libres (free shear layer ). Le
paragraphe suivant décrit l’évolution des régimes d’écoulements sur un cylindre en fonction
du nombre de Reynolds, afin de situer les cas tests étudiés. La Figure 2.5 de droite donne
une représentation schématique des différents régimes.
Ecoulement laminaire Pour un nombre de Reynolds inférieur à 4-5, le fluide est to-
talement attaché à la paroi (Figure 2.5 (b)). Jusqu’à des valeurs inférieures à 30-48, on
observe des recirculations symétriques laminaires dans le sillage proche (Figure 2.5 (c)).
L’allongement de ces recirculations provoque le démarrage d’une instationnarité (passage
de Figure 2.5 (c) à (d)). Les zones de cisaillement sont le siège d’une instabilité dites de
Kelvin-Helmoltz, due à un gradient de vitesse qui va engendrer une vitesse transverse à
l’écoulement. On assiste à une oscillation sinusoïdale des zones de cisaillement qui provoque
des détachements réguliers en opposition de phase de chaque côté du sillage. Les tourbillons
laminaires créés sont connus sous le nom d’allées de Bénard-Von Kármán.
Transition dans le sillage Pour la plage 200 < Re < 400, les tourbillons détachés
deviennent progressivement turbulents. Cette transition a lieu d’abord au loin, puis tend
à se rapprocher et à intervenir lors de la formation du tourbillon.
Transition dans les couches limites A des Reynolds plus élevés, la transition a lieu
dans la couche limite. La couche limite peut alors supporter plus longtemps le gradient
adverse de pression, ce qui décale le point de décollement vers l’aval. La pression exercée
sur l’arrière du cylindre est accrue, ce qui réduit drastiquement la poussée dans la direction
longitudinale (d’un facteur 4 environ). On retrouve ce phénomène sous le nom de drag crisis
dans la littérature (fig 2.5 (e)). C’est pour cette raison que les balles de golf, par exemple,
on des surfaces à aspérités pour se trouver dans ce régime et allonger les trajectoires. On
distingue dans ce régime une plage ’pré-critique’ et une plage ’post-critique’ qui encadrent
le phénomène de drag crisis
Régime pleinement turbulent Ce régime est atteint lorsque l’écoulement est turbulent
dans toute les régions entourant le cylindre. Le nombre de Reynolds pour lequel ce type
d’écoulement démarre est encore inconnu.
En terme de force, l’écoulement exerce sur le corps des forces statiques (Re < 30) ou
dynamiques. Il est habituel de les décomposer en une composante longitudinale et transver-
sale. On les nomme respectivement traînée (drag) et portance (lift). Il est usuel d’adimen-
sionnaliser ces forces par la force due à la pression dynamique sur une surface équivalente.
Selon la notation consacrée, on utilisera donc les coefficients de traînée CD et de portance
CL . On pose de même le nombre adimensionnel de Strouhal qui rend compte de la fréquence
des détachements.
1
F⊥ = ρSV 2 CL (2.32)
2
1
Fk = ρSV 2 CD (2.33)
2
fD
St = (2.34)
V
où :
• Fk et F⊥ sont les forces appliquées sur le cylindre respectivement dans les directions
longitudinale (parallèle) et transversale (orthogonale) à l’écoulement (N)
Les observations et calculs reportés par Noack [33] et Stansby [34], montrent que les
effets tridimensionnels de l’écoulement apparaissent au-delà d’une valeur du nombre de
Reynolds autour de 200. Cette valeur est particulièrement étudiée dans la littérature pour
des cylindres fixes. De plus, dans le cadre de la présente étude, nous nous intéresserons
au cas d’un cylindre en rotation. Nous prendrons comme référence le travail de Mittal et
Kumar [35]. Ils étudient par simulation numérique un cylindre en rotation uniforme dans
un écoulement laminaire pour des valeurs du paramètre α allant de 0 à 5. Le paramètre α
est similaire au paramètre d’avance d’une pale d’hydrolienne, il est défini comme le rapport
de la vitesse en paroi sur la vitesse au loin, voir (2.35).
Dans le cas du cylindre fixe, la portance suit une oscillation sinusoïdale de moyenne
nulle. Leur travail montre l’augmentation du coefficient de portance CL moyen avec la
vitesse de rotation. Ce phénomène est bien connu sous le nom d’effet Magnus et illustré
par les effets coupés et liftés au tennis. De façon intuitive, lorsque le cylindre tourne, le
fluide est accéléré d’un côté du fluide et ralenti de l’autre. On obtient donc une différence
de pression entre la partie supérieure et inférieure du cylindre. Les oscillations disparaissent
à partir d’une valeur de α = 1.91, puis réapparaissent pour α compris entre 4.3 et 4.75.
L’accent est mis sur l’hystérésis de l’instabilité obtenue par initialisation des simulations
pour des valeurs de α pour lesquelles on observe des allées de Karman.
Ω∗R
α= (2.35)
V
La comparaison s’appuiera également sur les travaux de James McNaugthon de l’uni-
versité de Manchester [30] qui a réalisé pendant sa thèse la validation d’une méthode sliding
mesh avec Code_Saturne. Il montre notamment que sa méthode nécessite plusieurs ité-
rations sur la résolution de Navier-Stokes à chaque pas pour obtenir la convergence, ainsi
qu’un déplacement à l’interface petit devant la taille des mailles. Ces deux contraintes
accroissent le temps de calcul. Nous cherchons donc quelles sont les limitations de notre
propre méthode.
Dans l’étude qui suit, on se sert de la méthode classique comme référence de calcul.
L’objectif de la validation est de reproduire le cas non-couplé à l’aide du cas couplé. La
simulation numérique donnant nécessairement un écart avec la réalité, on peut ainsi évaluer
l’erreur imputable au couplage uniquement.
Chacun des deux maillages est choisi structuré et conforme. A cause de la rotation, le
maillage n’est plus conforme à l’interface pendant la simulation. La situation se retrouve
sur la Figure 2.6 (b). On désigne par :
2.6. CAS TEST 1 : CYLINDRE 2D EN ÉCOULEMENT LAMINAIRE 59
Le couplage à l’interface est réalisé de telle sorte que le flux sortant d’une cellule A du
maillage Ωi soit transmis intégralement au maillage Ωj par la cellule B contenant le projeté
orthogonal du centre de gravité de A. Il convient donc que la taille des mailles à l’interface
soit de taille similaire. On choisit d’imposer le même nombre de segments de chaque côté.
L’étude de convergence en maillage utilise comme critère les coefficients de portance et
de traînée. Dans cette étude, on regarde la valeur moyenne du coefficient de traînée CD et
la moyenne quadratique du coefficient de portance CL sur un cylindre fixe. Les diffé-
Rext
rents paramètres sont reportés dans le tableau 2.1. Le rapport doit être suffisamment
Rc
important pour que le calcul dans un domaine fini reproduise le comportement en milieu
infini. Dans le cas contraire, on observe un effet de blocage qui accélère l’écoulement au
niveau du cylindre. La valeur du rapport varie entre les maillages A, B et C. La variation
des coefficients entre B et C est inférieure à 1%, on utilise donc le maillage B.
Les grandeurs utilisées ont pour valeurs :
• Rc = 0.025m
• Ri = 0.1m
60 CHAPITRE 2. LE MAILLAGE ROTATIF DANS CODE_SATURNE
• Rext = 2m
• U0 = 0.004m.s−1
• ρ = 1000kg.m−3
• µ = 10−3 kg.m−1 .s−1
Entre les maillages I à IV, on effectue un raffinement progressif de la taille des mailles,
et une convergence acceptable est observée à partir de BII. Le pas de temps de chaque
simulation est adapté afin de conserver un nombre de Courant maximum de l’ordre de
1. Le nombre de courant est approximativement le nombre de mailles traversées par une
particule fluide en un pas de temps. Le solveur utilise un schéma implicite, les calculs
convergeront donc même sans cette condition, mais elle est nécessaire afin d’optimiser la
précision du résultat dans chaque cas. Dans ces simulations, les restrictions sur le nombre
de maille à l’interface, imposées par la méthode de couplage, n’apparaissent pas. Lors du
passage au cylindre rotatif, on remarque que le calcul est convergé à partir de 180 mailles
sur la circonférence. Sauf mention particulière, on utilise le maillage B IV pour la suite de
l’étude.
0.8
0.6
0.4
0.2
−0.2
−0.4
−0.6
−0.8
100 120 140 160 180 200 220 240 260 280 300
Figure 2.7 – Coefficient de portance en fonction du temps pour les méthodes recollé (bleu)
et couplée (rouge)
Dans le cas couplé, sans rotation, on se retrouve avec une interface conforme. On élimine
dès lors le problème de l’absence d’interpolation pour le passage du flux de masse entre
cellules décalées. Pour vérifier cette propriété, on calcule la différence de flux de chaque côté
de l’interface. On impose de chaque côté la même condition de Dirichlet pour la vitesse.
Pour le calcul immobile, la propriété est vérifiée à l’ordre de la précision machine.
La Figure 2.7 compare les courbes de portance pour un calcul recollé et un calcul
utilisant la méthode de couplage. Le résultat du couplage est légèrement en avance sur la
méthode classique au démarrage de l’instabilité, mais en régime établi, amplitude, valeur
moyenne et fréquence sont identiques.
On peut conclure que le couplage fonctionne bien pour deux maillages immobiles et
conformes dans le cas d’une simulation laminaire.
62 CHAPITRE 2. LE MAILLAGE ROTATIF DANS CODE_SATURNE
∆t · Ω · nc
γ= (2.36)
2π
Table 2.2 – Erreur sur le flux de masse à l’interface. Les cases comportant * correspondent
aux calculs non-convergés
2.6. CAS TEST 1 : CYLINDRE 2D EN ÉCOULEMENT LAMINAIRE 63
Re= 200 alpha= 1.0 dt= 0.2s NTERUP= 1 ORDRE- 1 Re= 200 alpha= 1.0 dt= 0.2s NTERUP= 1 ORDRE- 1
8e- 04
defaut de m asse au passage de l interface % de m asse perdue entre inlet/ outlet
6e- 04
6e- 04
4e- 04
4e- 04
difference de flux adim ensionnelle
2e- 04
2e- 04
% de inlet
0e+ 00
0e+ 00
- 2e- 04
- 2e- 04
- 4e- 04 - 4e- 04
- 6e- 04
10000 10200 10400 10600 10800 11000 11200 11400 11600 11800 12000 10000 10500 11000 11500 12000
Iteration Tem ps
Figure 2.8 – Différence de flux au passage des conditions aux limites de couplage (a) et
masse perdue entre l’inlet et l’outlet. α = 1.0, au premier ordre en temps, ∆t = 0.2s
la conservation du débit entre l’entrée et la sortie est tracée en fonction du pas de temps
par la figure 2.8.b). Elle montre que les oscillations basse fréquence sont répercutées sur
l’ensemble mais que l’erreur locale sur le flux est compensée. Typiquement, l’erreur relative
est de 10−4 pour une précision sur la résolution des grandeurs de 10−6
Les simulations ont été réalisées avec des maillages comportant des nombres de mailles
P compris entre 120 et 480. On retrouve dans le tableau 2.2 les valeurs de l’écart type du
pourcentage de masse perdue à l’interface amont. Les calculs avec les maillages BII et BV
(P=N=480) divergent entre 8000 et 10000 itérations environ. La valeur reportée pour ces
cas ne tient donc pas compte de la partie divergente. Dans la partie divergente, l’erreur sur
le flux prend des valeurs jusqu’à cent fois supérieures au régime normal. Par contre pour
la masse totale perdue, les valeurs n’explosent pas. La masse perdue augmente jusqu’à
divergence du calcul.
L’erreur diminue avec l’augmentation du nombre de mailles. Dans le cas des maillages
BIII et BIV, le tableau 2 montre que la combinaison Rotor -Stator disposant respectivement
de 240 et 180 mailles sur la circonférence procure une erreur quasi constante, indépendante
du temps. Il est même étonnant de constater que ces valeurs sont toujours valables pour
des déplacements de l’ordre de 2 mailles par pas de temps (ie γ = 2). Le choix de ne
pas implémenter d’interpolation entre les faces frontières pour un meilleur passage du flux
parait donc se justifier. On remarque une augmentation à l’ordre 2. Cette donnée tend à
infirmer l’hypothèse que l’erreur est due au passage du flux à l’interface. En effet, le calcul
converge plus facilement avec le second ordre et donne de meilleurs résultats sur les efforts.
−1.6 −1.6
sans interface ordre 1 sans interface ordre 1
NTERUP=1 NTERUP=1
−1.8 NTERUP=2 −1.8 NTERUP=2
NTERUP=3
−2 −2
−2.2 −2.2
−2.4 −2.4
−2.6 −2.6
−2.8 −2.8
−3 −3
−3.2 −3.2
−3.4 −3.4
−3.6 −3.6
0 10 20 30 40 50 60 70 80 0 10 20 30 40 50 60 70 80
(a) γ = 1, 2 (b) γ = 0, 15
Figure 2.9 – Comparaison des coefficients de portance obtenus avec différents nombres
d’itérations
Champ de pression
α=1
ɤ = 2.44
Figure 2.10 – Discontinuité du champ de pression à l’interface Rotor /Stator aval, pour
une valeur de γ = 2.44
Il conviendra donc de ne pas dépasser cette valeur pour les calculs en couplage. Elle est
valable en calcul laminaire. On vérifiera dans un deuxième cas test ce qui advient en régime
turbulent.
−1
Methode recollee Code_Saturne − Ordre2
−1.2 Sliding − Code_Saturne − Ordre1
Sliding Manchester − Ordre2
−1.4
Mittal reference
−1.6 Sliding − Code_Saturne − Ordre2
−1.8
−2
−2.2
−2.4
−2.6
−2.8
−3
−3.2
−3.4
20 25 30 35 40 45 50 55 60 65 70 75 80
Figure 2.11 – Coefficient de portance en fonction du temps pour les différentes méthodes
de simulation, α = 1.0
référence dans le cadre de nos simulations. La géométrie utilisée reste la même que dans le
cas de l’écoulement laminaire. La condition limite d’entrée est un flux de masse constant
dans la direction x. Le rayon extérieur du domaine vaut 40 fois le rayon du cylindre et 10
fois celui de l’interface de couplage.
Une étude de convergence est menée avec un maillage recollé pour le cylindre fixe.
Pour le maillage BII adopté, la taille de la première maille à la paroi est fixée à 0.1mm.
Le maillage du Rotor est structuré et orthogonal. La direction radiale est découpée en 87
segments dont la taille évolue en progression géométrique de facteur 1, 25. La partie Stator
possède 120 segments radiaux, avec une progression de facteur 1, 20 vers l’extérieur. Les
deux maillages comptent 240 segments orthoradiaux.
2.7.2 Résultats
2 1,2
Expérimental Aoki
Expérimental Aoki
1,8 PANS Elmiligui
PANS Elmiligui
1 kwSST recollé
kwSST recollé
1,6 V2f recollé
V2f recollé
1,4
0,8
1,2
Cd moyen
Cl moyen
1 0,6
0,8
0,4
0,6
0,4
0,2
0,2
0 0
0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1 0 0,2 0,4 0,6 0,8 1 1,2
α α
(a) (b)
L’écoulement turbulent autour d’un cylindre a été simulé pour 4 vitesses de rotation.
Le pas de temps d’intégration temporelle est choisi à 0.25ms après une étude de conver-
gence avec le maillage recollé et une vitesse de rotation nulle. Cette valeur permet en outre
d’obtenir une valeur du nombre de Courant inférieur à 1. Le cas laminaire a montré l’im-
portance de la condition CFL au premier ordre. Pour un paramètre d’avance de α = 0.3,
on obtient un paramètre γ ≃ 0.15. On reste donc dans la limite déterminée au chapitre 2.6
d’une demi-maille de déplacement par pas de temps pour toutes les vitesses de rotation de
l’étude.
2.7. CAS TEST 2 : CYLINDRE TURBULENT 71
La Figure 2.12 montre les coefficients de portance et de traînée moyens pour les diffé-
rentes valeurs de paramètre d’avance. La courbe pleine provient de l’expérience de Aoki et
la courbe en pointillée représente les simulations de Elmiligui avec le modèle PANS. Pour
le coefficient de portance les résultats d’Elmiligui sont confondus avec les résultats expé-
rimentaux pour α ≤ 0.5. Après cette valeur la différence s’accentue fortement. La chute
locale de la portance et sa remontée ne sont pas décrits par le modèle PANS qui donne une
portance en constante augmentation. Pour le coefficient de traînée, le comportement est
identique, bien que la concordance des résultats pour α ≤ 0.5 soit moins nette que la pré-
cédente. La simulation avec les modèles RANS du Code_Saturne sous-estime le coefficient
de traînée CD d’un ordre de grandeur de 25% environ. En revanche, le modèle kω − SST
reproduit assez fidèlement le CL moyen alors que le modèle v 2 − f le sur- estime largement.
Comme pour le cas laminaire, on se sert des simulations en maillage recollé pour valider la
méthode de couplage.
Modèle kω − SST
La Figure 2.13 compile les résultats des simulations réalisées avec le modèle kω − SST .
Les courbes pointillées montrent les valeurs en recollé tandis que les points sont obtenus par
la méthode de couplage. Les écarts entre les deux approches sont en général importants,
sauf peut être pour la valeur moyenne du coefficient de portance (α ≤ 0.3). Pour α = 0.0
les valeurs moyennes des coefficients sont identiques avec les deux approches, par contre les
amplitudes sont complètement différentes. Les écarts concernant le coefficient de traînée
en sliding atteignent presque six fois la valeur de celui obtenu par le calcul recollé. Sur la
Figure 2.14 (a) on peut voir le champ de vitesse pour la simulation en recollé (gauche) et
en sliding (droite). On observe que le sillage sliding est moins allongé. Le régime des deux
écoulements semble être différent.
Pour α > 0.0, l’accord entre les deux calculs est acceptable pour les valeurs moyennes
du coefficient de portance (Figure 2.13 (b)) tant que α ≤ 0, 5.
Le coefficient de traînée ne suit pas la tendance à la baisse montrée expérimentalement
et que l’on retrouve en calcul recollé. La différence au niveau des amplitudes observées pour
le cylindre fixe se confirme en rotatif. Si les oscillations restent faibles en simulation recollée,
elles augmentent pour le calcul en couplage. Un facteur 10 est atteint sur l’amplitude du
coefficient de traînée pour α = 0.6. La figure 2.14 (b) présente le champ de vitesse dans
le sillage du cylindre pour la valeur α = 0.3. La différence de Reynolds apparent entre les
deux méthodes s’accroît avec la vitesse de rotation. Le sillage du calcul non-couplé prend
une allure qui repousse la transition turbulente des zones de cisaillement, alors que celui
du calcul couplé montre un enroulement complet de ces mêmes zones. De même, l’écart
de niveau pour la viscosité turbulente est encore plus marqué ; les plus hauts niveaux en
72 CHAPITRE 2. LE MAILLAGE ROTATIF DANS CODE_SATURNE
1,4 1,8
Sliding 1er ordre
1,6 Recollé 1er ordre
1,2
1,4
1
1,2
0,8
Cd moyen
Cl moyen
1
0,6 0,8
0,6
0,4
0,4
0,2 Sliding 1er ordre 0,2
Recollé 1er ordre
0 0
0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1
Alpha Alpha
0,35 2,5
Amplitude Cd
0,3
Amplitude Cl
2
0,25
1,5
0,2
0,15 1
0,1
0,5
0,05
0 0
0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1
Alpha Alpha
sliding sont 2 fois inférieurs à ceux de la méthode recollé (figure 2.14.c,d). Le calcul sliding
a donc tendance à prévoir un écoulement moins visqueux.
Modèle v 2 − f
1,2 2
Sliding 1er ordre Sliding 1er ordre
Recollé 1er ordre 1,8
Recollé 1er ordre
1
1,6
1,4
0,8
1,2
Amplitude Cd
Cl moyen
0,6 1
0,8
0,4
0,6
0,4
0,2
0,2
0 0
0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1
ALPHA ALPHA
a) CD moyen v 2 − f b) CL moyen v 2 − f
1 3
Sliding 1er ordre Sliding 1er ordre
0,9 Recollé 1er ordre Recollé 1er ordre
2,5
0,8
0,7
2
Amplitude Cd
0,6
Amplitude Cl
0,5 1,5
0,4
1
0,3
0,2
0,5
0,1
0 0
0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1
ALPHA ALPHA
c) CD Amplitude v 2 − f d) CL Amplitude v 2 − f
Les mêmes simulations sont effectuées avec le modèle RANS v 2 − f . La Figure 2.15
présente les résultats comparés des simulations en méthode de couplage et en maillage
recollé. On observe une meilleure cohérence des résultats. Pour α = 0.0, la comparaison
des deux méthodes pour un cylindre fixe montre une très bonne concordance. Les valeurs
moyennes des efforts sont identiques ; l’écart d’amplitude est de l’ordre de 10% pour le
coefficient de portance. La visualisation des champs de vitesse est tout aussi éloquente
2.7. CAS TEST 2 : CYLINDRE TURBULENT 75
(Figure 2.16 (a)) et illustre l’identité des écoulements. Pour α > 0.0, les valeurs α = 0.3 et
α = 0.6 montrent un accord très acceptable entre les deux méthodes, en particulier pour
la valeur moyenne du coefficient de portance. On obtient moins de 10% d’erreur sur le CD
moyen et 5% sur le CL moyen. Comme déjà mentionné à propos du modèle kω − SST ,
la valeur des amplitudes est bien plus grande en configuration sliding, qu’en configuration
recollé. Les Figures 2.16 (b) et (d) illustre respectivement le champ de vitesse et de viscosité
turbulente pour α = 0.3. On remarque de nouveau Reynolds apparent plus élevé en sliding
ainsi qu’une perte de 50% de la viscosité turbulente, ce qui explique l’écart au niveau des
amplitudes d’oscillation. Il semble de plus que la viscosité turbulente ’adhère’ à l’interface
et est entraînée dans le sens de la rotation. Pour des valeurs de α plus élevées, on obtient
une augmentation des différences. Le CD moyen augmente, contrairement à la tendance
observée. Le CL moyen reste proche de la méthode recollé, mais son amplitude devient 5
fois plus élevée. On remarque que l’amplitude des oscillations de traînée augmente dans les
mêmes proportions sur cette dernière simulation, tout en restant proche du recollé.
76 CHAPITRE 2. LE MAILLAGE ROTATIF DANS CODE_SATURNE
Modèle k − ǫ
Les problèmes rencontrés avec les modèles bas Reynolds, et spécialement le modèle
kω − SST amènent à étudier le comportement du modèle k − ǫ. Il s’agit en quelque
sorte de simplifier le problème et de minimiser les sources d’erreur. L’implémentation du
modèle k − ǫ est assez standard et ainsi plus fiable. On utilise un maillage similaire aux cas
précédents, quoi que modifié en proche paroi ; le seul paramètre à varier est la taille de la
première maille à la paroi. Une taille de 0.5mm permet d’obtenir une valeur y+ maximum
comprise entre 25 et 40.
Les Figures 2.17 a-d) illustrent, de façon similaire aux autres modèles, les champs de
vitesse et de viscosité turbulente. Comme pour le modèle v 2 − f , dans le cas du cylindre
fixe, on observe un très bon accord entre les méthodes couplé et recollé. Pour un paramètre
d’avance non nul (α = 0.6 sur les figures 2.17 b et d ), la chute de viscosité turbulente est
encore constatée. Par rapport au précédents modèles, la viscosité turbulente est plus forte
et fige les instationnarités comme bien visible en 2.17.d.
2.8.1 Maillages
a) b)
Code_Saturne Comme cela a été dit, cette étude fait suite à une première étape de
validation de la méthode de maillage glissant de Code_Saturne pour des écoulements
laminaires et turbulents autour d’un cylindre. Elle a mis en évidence le dysfonctionnement
2.8. SIMULATION DE L’HYDROLIENNE ACHARD10 EN 2D 81
de la technique implémentée dans le code avec le modèle kω − SST , mais une concordance
acceptable du calcul est obtenue en modèle v 2 −f par rapport aux résultats attendus. Nous
utilisons donc ce dernier modèle pour effectuer la simulation de l’hydrolienne en 2D. Un
calcul à l’optimum avec le modèle kω − SST , permettra d’estimer la différence entre les
deux modèles.
Le modèle bas-Reynolds v 2 − f propose une modélisation de la sous-couche visqueuse.
La première maille en paroi doit donc se situer en-dessous d’une valeur de y + = 2, 3. Le
maillage kω − SST élaboré par Ervin Amet est donc également adapté à un calcul en
v2 − f .
2.8.3 Résultats
La Figure 2.19 présente les coefficients de puissance (CP ) moyens obtenus par calcul
2D avec Code_Saturne et Ansys-Fluent, ainsi que les résultats expérimentaux. Les calculs
2D surestiment de façon importante le CP pour les paramètres d’avance supérieurs à 1.
Dans la partie à grands λ, les trois courbes sont parallèles, ce qui montre une assez bonne
reproduction des efforts par le calcul lorsque les décollements disparaissent. On conserve
un écart dû aux effets 3D. Néanmoins, le calcul avec Code_Saturne surestime moins la
valeur expérimentale. Il a été montré dans [41] que ces effets 3D sont dus essentiellement
aux tourbillons créés par le raccord bras-pale et aux vortex d’extrémité de pale. Les pertes
par frottement sur les bras sont négligeables.
A petit λ, le CP donné par Ansys-Fluent est inférieur à la valeur obtenue dans la
veine d’essai, alors que Code_Saturne affiche une valeur similaire (Figure 2.19). On notera
l’influence de deux phénomènes relatifs au cas 3D qui ne sont pas reproduits par le calcul
3D :
◦ Un effet positif sur le CP existe par le passage en 2D car on ne prend pas en compte
les tourbillons marginaux de bout de pale ainsi que la traînée exercée sur les bras de
raccord.
0.7
0.65
0.6
0.55
0.5
0.45
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0
1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3
la part du coefficient de puissance fournie par la pale numéro 1. Cette grandeur n’est pas
mesurable expérimentalement par notre technique. Les Figures 2.20, 2.24 et 2.25 présentent
les courbes de CP instantanées calculées ainsi que les valeurs expérimentales (sauf pour
2.20.b) minimales et maximales obtenues au moyen de l’évaluation des incertitudes de
mesure.
De manière générale, on remarque que les courbes sont de plus en plus régulières avec
l’augmentation du paramètre d’avance. De même, la valeur moyenne du CP est convergée
au bout de 3 tours pour λ =3, alors que cette convergence n’est pas atteinte à moins de
5% pour λ =1. La raison principale est la présence des forts décollements qui régissent
l’écoulement et leur déclenchement aléatoire sur une plage de θ. La forme des courbes
obtenues par les mesures et l’expérience pour λ =2 est globalement la même alors qu’elle
est très différente pour le régime sub-optimal. L’écoulement calculé est probablement assez
différent de l’écoulement réel.
Les Figures 2.21,2.22 et 2.21 présentent les champs de vorticité suivant la verticale (la
composante perpendiculaire au plan 2D) mesurés et calculés autour des pales de l’hydro-
lienne, respectivement pour λ =2 et 1. Si l’écoulement est suffisamment bi-dimensionnel,
cette composante est la seule non-nulle. La thèse de Jonathan Bossard [5] montre que cette
hypothèse est bien vérifiée dans le plan mi-pale pour des vitesses spécifiques supérieures à
λ =1,75.
L’analyse détaillée des résultats expérimentaux ainsi que la comparaison avec les si-
mulations Ansys-Fluent sont présentées dans le manuscrit de thèse de Jonathan Bossard
2.8. SIMULATION DE L’HYDROLIENNE ACHARD10 EN 2D 83
1.2 0.8
1.1
0.7
1
0.9 0.6
0.8
0.5
0.7
0.6 0.4
0.5
0.4 0.3
0.3
0.2
0.2
0.1 0.1
0
0
−0.1
−0.2 −0.1
0 50 100 150 200 250 300 350 400 0 50 100 150 200 250 300 350 400
(Chap.4, p.144-157). On ne donnera ici qu’un résumé de cette analyse ainsi que les conclu-
sions majeures qui s’en dégagent. Dans le présent travail, on s’intéressera particulièrement
aux différences qu’apporte la simulation avec Code_Saturne et le modèle v 2 − f .
λ =2
Les mesures PIV nous permettent les observations suivantes : pour le premier tiers
(θ < 120◦ ) de la rotation, l’écoulement reste attaché au profil. C’est une mise en évidence
du retard au décollement dû aux effets dynamiques. La position θ = 120◦ correspond à un
angle d’attaque de α = 30◦ alors que l’angle de décrochage statique se situe aux alentours de
α = 10◦ [43]. Notons cependant que l’angle d’incidence sur la pale de la turbine est calculé
avec une vitesse axiale égale à U∞ alors qu’à l’approche de la turbine, la vitesse axiale
diminue sensiblement (d’1/3 au maximum suivant le théorème de Betz) ce qui a pour effet
de réduire l’angle effectif d’attaque sur la pale. A partir de θ = 140◦ (α = 27◦ ), l’écoulement
subit un décollement caractérisé par une recirculation au bord de fuite (vorticité négative,
bleue). A θ = 180◦ , les deux vortex grossissent et le tourbillon de décrochage dynamique se
détache. A partir de θ = 200◦ , les 2 tourbillons se détachent dans le sillage. Ils interagissent
avec les pales suivantes, à la position θ = 220◦ environ. A partir de θ = 240◦ , l’écoulement
redevient attaché. On n’observe pas d’autres détachement durant la fin du tour.
La simulation Ansys-Fluent montre un décrochage dynamique initié à partir de θ =
140◦ , donc en retard de 20◦ sur l’expérience. Le phénomène de décrochage se poursuit
jusqu’à θ = 220◦ . A cette position les deux tourbillons contrarotatifs sont lâchés dans le
84 CHAPITRE 2. LE MAILLAGE ROTATIF DANS CODE_SATURNE
Figure 2.21 – Champs de vorticité autour d’une pale à λ =2. Comparaison entre les
mesures PIV de Jonathan Bossard [5], les simulations numériques de Jéronimo Zanette [3],
les simulations Code_Saturne
2.8. SIMULATION DE L’HYDROLIENNE ACHARD10 EN 2D 85
Figure 2.22 – Champs de vorticité autour d’une pale à λ =2. Comparaison entre les
mesures PIV de Jonathan Bossard [5], les simulations numériques de Jéronimo Zanette [3],
les simulations Code_Saturne
86 CHAPITRE 2. LE MAILLAGE ROTATIF DANS CODE_SATURNE
Figure 2.23 – Champs de vitesse relative autour de la pale à différentes positions pour
λ =2. Les vecteurs vitesses sont représentés avec une norme constante, et superposés au
champ de vorticité.
2.8. SIMULATION DE L’HYDROLIENNE ACHARD10 EN 2D 87
sillage de la pale. L’hypothèse énoncée pour ce retard est l’état turbulent de la sous-couche
visqueuse dans la simulation, que l’on ne trouve pas dans l’expérience.
Les résultats de la simulation v 2 − f avec Code_Saturne, montrent une apparition
anticipée du gradient de pression adverse θ = 110◦ , qui se transforme en recirculation,
comme le montre la Figure 2.23. On observe une migration lente du point de décollement
le long de l’intrados. Le lâcher complet des tourbillons intervient à partir de θ = 200◦ .
On observe ce phénomène sur le CP , avec une inflexion de la pente montante du pic plus
précoce que sur celle du résultat Ansys-Fluent. Ensuite la chute de puissance est aussi plus
lente.
On observe que la forme des structures lâchées en v 2 − f est plus proche de la visua-
lisation PIV. De plus, le point de détachement en v 2 − f se situe à environ au milieu de
l’intrados et dans le dernier quart pour le modèle kω − SST . Ce point se situe au milieu
de la corde expérimentalement comme on peut le voir à θ = 180◦ .
λ =1
Pour le paramètre d’avance le plus faible (λ =1), la formation du tourbillon de décro-
chage dynamique commence vers θ = 80◦ . De plus, ce tourbillon reste au bord d’attaque
jusqu’à θ = 180◦ . En effet, la vitesse de rotation étant égale à celle de l’écoulement, la
pale ne peut dépasser le tourbillon. Au delà de cette position, le tourbillon est chassé à
l’extrados, et le très fort angle d’incidence induit des décollements au bord de fuite et au
bord d’attaque, qui sont lâchés respectivement à θ = 230◦ et θ = 240◦ .
La dynamique est sensiblement bien reproduite par les simulations 2D. On observe pour-
tant des différences marquées. Entre θ = 60◦ et 80◦ , la double recirculation présente dans
l’expérience et la simulation Ansys-Fluent ne se retrouve pas sur le calcul Code_Saturne.
Ensuite, à partir de θ = 100◦ , les deux modèles montrent l’apparition d’une recirculation
dans le sillage qui remonte jusqu’au bord de fuite et à l’intrados. Ce phénomène qui n’est
pas présent lors des mesures PIV persiste jusqu’à θ = 190◦ avec une forte dissipation pour
Code_Saturne et continue jusqu’à θ = 200◦ pour Ansys-Fluent. Le tourbillon de bord
d’attaque reste attaché jusqu’à θ = 260◦ pour les deux modèles.
88 CHAPITRE 2. LE MAILLAGE ROTATIF DANS CODE_SATURNE
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05
−0.05
−0.1
−0.15
0 50 100 150 200 250 300 350 400
Figure 2.24 – Coefficients de puissance total sur un tour pour différents modèles, λ =1
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
−0.1
−0.2
−0.3
−0.4
0 50 100 150 200 250 300 350 400
Figure 2.25 – Coefficients de puissance total sur un tour pour différents modèles, λ =3
2.8. SIMULATION DE L’HYDROLIENNE ACHARD10 EN 2D 89
Figure 2.26 – Champs de vorticité autour d’une pale à λ =1. Comparaison entre les
mesures PIV de Jonathan Bossard [5], les simulations numériques de Jéronimo Zanette [3],
et les simulations Code_Saturne
90 CHAPITRE 2. LE MAILLAGE ROTATIF DANS CODE_SATURNE
Points clés
• Le code de CFD open source Code_Saturne est utilisé comme outil de simulation
pour les travaux de cette thèse.
• Le cas test d’un cylindre en rotation dans un écoulement laminaire présente une bonne
stabilité, ainsi que des résultats sur les efforts identiques aux données expérimentales,
et à une méthode numérique similaire.
• Plusieurs modèles à viscosité turbulente à deux équations sont testés pour la simu-
lation d’un cylindre en rotation dans un écoulement de Reynolds Re = 60k. Dans
tous les cas, la rotation avec interface de couplage sous estime la production turbu-
lente et abaisse la viscosité apparente en aval de l’obstacle. Un résultat similaire est
constaté avec le code commercial Fluent au premier ordre en temps. Ce problème
semble disparaître avec un calcul effectué au second ordre en temps.
• Une incohérence entre les calculs avec et sans interface de couplage est constatée avec
le modèle kω − SST , même dans le cas sans rotation. Ce constat amène à utiliser le
modèle v 2 − f pour la suite des simulations.
Bibliographie
[1] G. Brochier, P. Fraunie, C. Beguier, and I. Paraschivoiu, “Water channel experiments
of dynamic stall on darrieus wind turbine blades,” J Propuls, 1986.
[2] E. Amet, “Simulation numérique d’une hydrolienne à axe vertical de type darrieus,”
2008.
[6] A. Mentxaca, “Analyse numérique des hydroliennes à axe vertical munies d’un caré-
nage,” 2011.
[7] V. Aumelas, Modélisation des hydroliennes à axe vertical libres ou carénées : développe-
ment d’un moyen expérimental et d’un moyen numérique pour l’étude de la cavitation.
PhD thesis, Grenoble-INP, 2011.
[8] F. Archambeau, N. Mehitoua, and M. Sakiz, “Code saturne : A finite volume code for
turbulent flows,” Int. J. Finite Volumes, 2004.
[9] “Code_saturne : Edf’s general purpose cfd software, user meeting,” 2013.
[10] S. Rolfo, J. Uribe, and D. Laurence, “Les and hybrid rans/les of turbulent flow in fuel
rod bundle arranged with a triangular array,” in Direct and Large-Eddy Simulation
VII, pp. 409–414, Springer, 2010.
[16] P. Durbin, “Near-wall turbulence closure modeling without damping functions,” Theo-
retical and Computational Fluid Dynamics, vol. 3, pp. 1–13, 1991.
[19] J. Benek, “A flexible grid embedding technique with application to the euler equa-
tions,” in 6th Computational Fluid Dynamics Conference Danvers, 1983.
[20] J. Steger, “The chimera method of flow simulation,” in Workshop on applied CFD,
Univ of Tennessee Space Institute, 1991.
[22] R. Steijl and G. Barakos, “Sliding mesh algorithm for cfd analysis of helicopter rotor–
fuselage aerodynamics,” International journal for numerical methods in fluids, vol. 58,
no. 5, pp. 527–549, 2008.
[23] E. L. Blades and D. L. Marcum, “A sliding interface method for unsteady unstructured
flow simulations,” International journal for numerical methods in fluids, vol. 53, no. 3,
pp. 507–529, 2007.
[24] O. Petit, M. Page, M. Beaudoin, and H. Nilsson, “The ercoftac centrifugal pump open-
foam case-study,” in 3rd IAHR International Meeting of the Workgroup of Cavitation
and Dynamic Problems in Hydraulic Machinery and Systems, pp. 523–532, 2009.
[25] C. S. Peskin, “The immersed boundary method,” Acta numerica, vol. 11, 2002.
BIBLIOGRAPHIE 93
[26] A. Jendoubi, D. Yakoubi, A. Fortin, and C. Tibirna, “An immersed boundary method
for fluid flows around rigid objects,” Int. J. Numer. Meth. Engng, vol. 1, p. 3, 2012.
[28] www.code-saturne.org www.code saturne.org, “Code saturne 2.2.0 theory and pro-
grammer’s guide,” tech. rep., EDF R&D, 2011.
[29] J. Mc Naughton, Turbulence modelling in the near field of an axial flow tidal turbine
using Code_Saturne. PhD thesis, University of Manchester, 2013.
[31] M. M. Zdravkovich, Flow around circular cylinders, vol. 1. Oxford University Press,
1997.
[32] A. Roshko, “Experiments on the flow past a circular cylinder at very high reynolds
number,” Journal of Fluid Mechanics, no. 10, pp. 345–356, 1961.
[34] P. Stansby and R. Rainey, “A cfd study of the dynamic response of a rotating cylinder
in a current,” Journal of Fluids and Structures, vol. 15, no. 3–4, pp. 513–521, 2001.
[35] S. Mittal and B. Kumar, “Flow past a rotating cylinder,” Journal of Fluid Mechanics,
vol. 476, no. 4, pp. 303–334, 2003.
[37] P. T. Tokumaru and P. E. Dimotakis, “The lift of a cylinder executing rotary motions
in a uniform flow,” Journal of Fluid Mechanics, vol. 255, pp. 1–10, 1993.
[38] T. Ito and K. Aoki, “Flow characteristics around a rotating cylinder,” 9th International
Symposium on flow visualization, vol. 26, pp. 29–34, 2001.
94 CHAPITRE 2. LE MAILLAGE ROTATIF DANS CODE_SATURNE
[39] A. Elmiligui et al., “Numerical study of flow past a circular cylinder using rans, hybrid
urans/les and pans formulations„” AIAA Paper, 2004.
[40] A. C. Benim et al., “Rans predictions of turbulent flow past a circular cylinder over
the critical regime,” in In Proceedings of the 5th IASME / WSEAS, vol. 5, p. 232, aug
2007.
[41] T. Maitre, E. Amet, and C. Pellone, “Modeling of the flow in a darrieus water turbine :
wall grid refinement analysis and comparison with experiments,” Jour of Renewable
Energy, 2012.
[42] T. Maitre, E. Amet, and C. Pellone, “Modeling of the flow in a darrieus water turbine :
wall grid refinement analysis and comparison with experiments,” Jour of Renewable
Energy, accepted, 2012.
A l’heure de la mise en eau de nombreux prototypes sur les sites d’essais, et à l’aube
du développement industriel des hydroliennes, la compréhension de l’interaction entre ma-
chines placées dans un parc est un sujet d’actualité. D’un point de vue expérimental, les
modèles réduits de turbine ont permis depuis les années 90 des études très complètes des
machines seules. En ce qui concerne les parcs, les moyens d’essais en laboratoire d’une
ampleur suffisante sont peu nombreux [1, 2, 3]. Ils restent limités à quelques machines
de géométrie réaliste, ou utilisent une représentation des machines par des disques solides
poreux. Cette dernière solution reproduit avec fidélité l’écoulement derrière une turbine à
partir de 10 diamètres à l’aval environ [4] . L’équivalence de cette représentation n’est pas
prouvée pour les machines à flux transverse. Le calcul complet de l’écoulement hydrody-
namique autour des turbines est bien maîtrisé pour les deux types de géométrie.
Des simulations à l’aide des modèles de turbulence RANS, puis LES apportent une
connaissance accrue de l’interaction d’une machine seule avec l’écoulement. Cependant, le
temps de calcul que demanderait la simulation d’un parc en décrivant la géométrie complète
de toutes les machines reste prohibitif. De plus, on peut considérer que les différentes
échelles en jeu (structures tourbillonnaires, pales, hauteurs d’eau, bassin côtier) ne sont
pas accessibles dans un même modèle. Il apparaît nécessaire d’avoir des modèles simples
permettant de calculer et d’optimiser un parc d’hydroliennes.
Ce chapitre présente une revue des modèles existants pour la simulation de parc d’hy-
droliennes en mettant l’accent sur les solutions plus spécifiquement adaptées pour des
machines à flux transverse. Une bonne compréhension de ces modèles et de leur champ
d’application nous permet d’introduire, dans un premier temps, les apports et les atouts
de l’approche proposée dans ce chapitre. Dans un deuxième temps, la méthodologie de
construction du modèle de turbine est présentée, en particulier ses différents degrés de
liberté et la méthodologie utilisée pour les fixer. Les résultats en termes de puissance ré-
95
96 CHAPITRE 3. MÉTHODE D’ÉQUIVALENCE EMPIRIQUE
cupérée sur une machine seule, puis un couple de machines sont validés par rapport aux
simulations effectuées en géométrie complète.
T 1
= ρu2 (1 − α42 ) (3.1)
A 2
On établit l’équation de conservation de la quantité de mouvement à l’intérieur du volume
de contrôle. On peut montrer, en considérant un volume de contrôle suffisamment grand
que l’intégrale des pressions sur le contour est nulle. La variation d’énergie entre l’entrée
et la sortie est due à la force T :
Ui3
Ui2 Ui4
Ui1
A1
A
A4
1 2 4
3
Figure 3.1 – Représentation de la veine de courant pour une turbine libre en fluide parfait
En égalisant les deux expressions on obtient la relation entre les coefficients de vitesse :
α2 = 1+α2
4
. La vitesse dans le disque est la moyenne des vitesse amont et aval. La puissance
extraite est donc égale à P = α2 uT que l’on peut exprimer en fonction de α2 :
1 1
P = ρAu3 α2 (1 − α42 ) = ρu3 4α22 (1 − α2 ) (3.3)
2 2
Le coefficient de puissance adimensionnel est égal à :
P
CP = 1 = 4α22 (1 − α2 ) (3.4)
2
ρu3 A
En différenciant cette expression, on trouve que la puissance maximum récupérable est
égale à 16/27e de la puissance cinétique du fluide à vitesse u traversant une surface A. Ce
maximum est atteint pour un facteur α2 = 2/3. On remarquera que le coefficient de trainée
CT correspondant est égal à 8/9e .
Les conclusions fondamentales de cette analyse sont les suivantes :
1. Pour un système hydro-cinétique dans un écoulement libre, on ne peut récupérer plus
de 16/27e de la puissance incidente. Elle est généralement appelée limite de Betz. Ce
facteur est à ne pas confondre avec le rendement du système qui caractérise lui sa
capacité à transformer cette énergie récupérable en énergie mécanique profitable.
2. Pour un fonctionnement de la turbine à l’optimum, la vitesse dans la machine et
reliée théoriquement à la vitesse de référence amont.
Il est nécessaire de s’attarder sur les limites de cette théorie. Premièrement, la génération
de couple impose une réaction tangentielle sur la vitesse qui induit une rotation de la
98 CHAPITRE 3. MÉTHODE D’ÉQUIVALENCE EMPIRIQUE
veine à l’aval. Elle produit une diminution de la puissance récupérable. A puissance égale,
la vitesse de rotation diminue le couple nécessaire exercé sur les pales. La limite réelle
de puissance récupérable dépend du couple. Cette limite décroit vers les basses vitesses
spécifiques. Lors de la construction des machines axiales, on cherche à augmenter le plus
possible le point de fonctionnement optimum afin de diminuer les efforts de pression. Enfin,
un phénomène de mélange entre l’intérieur et l’extérieur a en réalité lieu, réalimentant la
vitesse du fluide dans le sillage, ce qui met en défaut l’application du théorème de Bernoulli
dans la partie aval. L’expérience a montré que les résultats précédents sont globalement
valides pour un coefficient d’induction α2 inférieur à 0, 6 environ c’est à dire en dessous de
l’optimum [8]. En réalité, le coefficient de trainée dépasse l’unité, ce qui entraînerait une
vitesse négative derrière la turbine dans la théorie de Betz. Ces mesures sont expliquées
par le mélange tourbillonnaire. Lorsque la différence u et u2 devient trop importante, une
zone de cisaillement apparaît ce qui va induire des instabilités favorisant la diffusion de
quantité de mouvement entre le fluide intérieur et extérieur.
Ue4 p
4
p1
Ui3
Ui2
Ui1 Ui4
U5
A1 p3
p2 A
p5
A4
1 2 4
3
(a) (b)
celle dans laquelle la limite stricte (par rapport à la vitesse dans le canal sans turbine) est
dépassée. Dans la zone grise, la limite de Betz par machine ne peut être atteinte.
Une autre famille des modèles analytiques est celle des modèles utilisant la Blade Ele-
ment Momentum Theory (BEMT). Cette méthode couple le calcul des forces sur la turbine
par la méthode des éléments de pale et la théorie du disque d’action, qui égalise la force ap-
pliquée sur les pales avec la quantité de mouvement perdue par le fluide. Cette technique
permet d’obtenir un calcul des performances globalement satisfaisant dans le cas d’une
machine seule en écoulement libre. Son utilisation pour les machines à flux transverse est
analysée en détail par Islam et. al. [16] mais également dans les travaux de Gretton [17]. La
dynamique complexe de l’écoulement a amené les chercheurs à élaborer des modèles évo-
lués. La veine de courant est découpée en plusieurs tubes, qui en suivent l’expansion [18].
L’extraction d’énergie est répartie en deux zones d’action, amont et aval. Enfin les pertes
par détachement tourbillonnaires et par tourbillons de bout de pales sont prises en compte
par des méthodes de vortex. Si ces méthodes permettent de retrouver les performances
mesurées, elles restent très sensibles aux données d’entrée sur les profils des pales et aux
effets instationnaires. De plus, elles ont des difficultés pour représenter des rotors à forte
solidité. Enfin, le fondement analytique de ces méthodes fait qu’elles ont par construction
les mêmes limites que les modèles idéaux présentés précédemment. Ainsi, elles ne sont pas
utilisées pour les calculs de performance en configurations de parc.
102 CHAPITRE 3. MÉTHODE D’ÉQUIVALENCE EMPIRIQUE
Pertes de charges dans une simulation d’écoulement visqueux Une seconde stra-
tégie qui est privilégiée actuellement pour les calculs de parc dans des environnements
réalistes est l’utilisation du couplage entre une simulation d’écoulement et des termes de
quantité de mouvement. Dans cette section, nous pouvons faire la distinction entre les mo-
dèles qui utilisent des valeurs de performance empiriques et ceux qui en font le calcul au
sein de la simulation.
Les modèles dits RANS+BE intègrent les calculs des efforts par la théorie des éléments
de pales, dans une simulation d’écoulement classique avec modèles de turbulence RANS. Ce
type de couplage est devenu presque classique pour les machines à flux axial [20, 21, 22, 23],
mais reste encore d’utilisation limitée pour celle à flux transverse. Le travail de Antheaume
et al., réalisé au LEGI, exploite cette méthode pour étudier une ligne de machines dans
un écoulement. Il utilise un maillage 3D pour l’écoulement global, avec une construction
spécifique dans l’environnement de la machine (voir Figure 3.4). Le trajet parcouru par
les pales est discrétisé sur une couronne. Pour chaque maille, on calcule la vitesse relative
3.1. LES PARCS ET LES MODÈLES EXISTANTS 103
du fluide vue par la pale et ainsi son l’angle d’incidence, à partir de la simulation RANS.
On peut en déduire la valeur de la force appliquée sur la pale, en moyenne sur un tour, à
partir des courbes de portance et de trainée empiriques pour le profil de pale choisi. Une
force équivalente est ensuite reproduite dans l’écoulement via un terme source/puits de
quantité de mouvement. La méthode donne de bons résultats sur le coefficient de puissance
à partir de la région de l’optimum. Les écarts s’expliquent par la prépondérance des effets
tourbillonnaires, et d’hystérésis non considérés dans le modèle.
Une option plus simple mais basée sur la théorie analytique est celle d’utiliser des
données empiriques pour les coefficients de force et de puissance. Dans ce cas, on utilise
souvent les données expérimentales corrigées du blocage de la section de test pour se
ramener au cas infini [24]. La force de trainée est dans ce cas également représentée par un
terme de perte de charge dans l’écoulement, que l’on calcule grâce au coefficient Ct,exp et
à la vitesse locale, U∞ la vitesse de l’écoulement libre, ou de l’écoulement amont dans un
canal. Par définition du coefficient adimensionnel, la force de trainée s’exprime en fonction
de la vitesse de l’écoulement libre non contraint par :
1 2
Ft = − ρAt Ct,exp U∞ (3.6)
2
Dans le cas d’un parc, la présence d’autres turbines crée un effet de confinement qui ne
permet pas d’avoir cette vitesse de référence. Les modèles utilisent généralement les résul-
tats de la théorie de Betz 1D, généralisée dans GC07, qui relient la vitesse de référence à
la vitesse locale via le facteur d’induction a. Le facteur d’induction est ensuite pris égal à
1/3 (optimum de puissance) ou dérivé de la formule :
Ulocal
U∞ = (3.7)
(1 − a)
1 p
a = (1 − 1 − Ct ) (3.8)
2
p
1 1 − 1 − Ct,exp
Ft = − ρAt (4 · p )|Ulocale |Ulocale (3.9)
2 1 + 1 − Ct,exp
Ces méthodes utilisent donc l’hypothèse forte que le facteur d’induction est constant
quelque soit la configuration. En d’autres termes, cela revient à considérer que l’effet de
blocage est directement équivalent à une augmentation de la vitesse U∞ , dans une situation
où la machine est supposée être déconfinée. Ce résultat théorique, valide pour un générateur
idéal dans un écoulement sans perte, sera discuté dans ce chapitre. Le même déroulement
analytique permet de connaître la puissance en fonction du coefficient de puissance ex-
opt
périmental CP,exp et la vitesse locale. Roc et al. [25] appliquent cette méthode couplée à
104 CHAPITRE 3. MÉTHODE D’ÉQUIVALENCE EMPIRIQUE
W Centre O1
λ1 , q1X , CP1
V0 Centre Om
λm , qmX , CPm
Centre OM
L λM , qMX , CPM
Figure 3.5 – Schéma de principe d’un calcul d’hydrolienne dans un domaine quelconque,
résumé des notations
la suite la construction des données empiriques et les tests sur leur validité. Le principal
paramètre dont s’affranchit le modèle est le Nombre de Reynolds de l’écoulement. Les me-
sures effectuées durant la thèse de V. Aumelas [28] ont montré qu’au delà d’une valeur
critique de Re, son influence sur les puissances et trainées relatives est faible. La méthode
dans son implémentation première détaillée ici, se base sur le calcul d’un écoulement bidi-
mensionnel afin de valider “l’équivalence” entre la représentation par pertes de charge et la
présence d’une vraie machine. Il n’y a pas de difficulté identifiée pour un passage en 3D. Il
s’agira seulement d’établir les courbes empiriques pour une machine en trois dimensions.
fluide. L’extraction d’énergie par la turbine, crée un saut de pression entre l’amont et
l’aval. Ce saut de pression équivaut à une force Fm exercée par le fluide sur la turbine, dont
les composantes sont désignées par Fm,X la trainée globale et Fm,Y la force transverse. Pm
désigne la puissance extraite par la turbine m. Cette force est représentée dans l’écoulement
par deux distributions spatiales termes sources de quantité de mouvement qmX et qmY . On
a donc la relation suivante entre les distributions qm et la force exercée sur la turbine :
ZZ
FmX,Y = −h qmX,Y dS (3.10)
S
La méthode utilise une approche stationnaire. Ainsi, les grandeurs décrites par la suite
représentent leur valeur moyennée sur un tour (ou 1/3 de tour). A partir de ces forces,
nous pouvons définir des coefficients adimensionnels locaux, c’est à dire prenant comme
référence le débit d’énergie cinétique traversant effectivement le disque. Il est à noter que
cette référence est plus faible que celle habituellement considérée, ie le débit d’énergie
cinétique traversée par l’écoulement libre. Ces coefficients, que l’on appellera locaux sont
∗
notés avec le symbole (∗). On construit les coefficients locaux de trainée CmX,Y :
∗ FmX,Y
CmX,Y = (3.11)
ρRhVm2
le coefficient de puissance local
Pm
CP∗m = (3.12)
ρRhVm3
La construction des données empiriques revient à établir la relation entre ces coefficients
locaux et le paramètre d’avance local λ∗ .
λ∗m = ωm R
Vm
(3.13)
Dans le disque de la turbine, le couple [qmX , qmY ] de distributions spatiales sont deux
fonctions choisies arbitrairement par similarité phénoménologique. L’écoulement est calculé
par une méthode de discrétisation volume fini dans Code_Saturne avec le modèle haut-
Reynolds k−ǫ. Une implémentation simple permet d’ajouter les termes sources à l’équation
de conservation de la quantité de mouvement :
∂ (ρ~v ) →− →
− →
−
+ ∇ · (ρ~v ⊗ ~v ) = − ∇p + ∇ · τ + −
q→
m (3.14)
∂t
Etablir le point de fonctionnement pour le parc de M machines nécessite de déterminer
les triplets [FmX , FmY , Vm ] pour chacune des turbines. On procède par une méthode itéra-
tive jusqu’à converger vers un état stable en moyenne. La procédure est résumée par les
points suivants :
108 CHAPITRE 3. MÉTHODE D’ÉQUIVALENCE EMPIRIQUE
n−1
1. Calcul de l’écoulement avec les valeurs de FmX,Y
2. Détermination de Vm et λ∗m
n
3. Détermination FmX,Y et CP m par utilisation des données empiriques au pas de temps
n et passage des valeurs pour l’étape n + 1
S1 S1
S S2 S S2
Figure 3.6 – Champ de vitesse et lignes de courant autour de l’hydrolienne pour deux
situations de confinement, ǫ = 2 à gauche, et ǫ = 5, 7, à droite
Données empiriques issues du calcul L’objectif du modèle est de faire un calcul des
efforts et des puissances pour des machines à l’intérieur d’une modélisation de l’écoulement.
Il doit nécessairement être simple et adaptable à d’autres machines après une étude simple
des performances, numérique ou expérimentale. Cette partie a pour objectif de présenter
les calculs effectués sur la machine complète dans différentes configurations et d’expliquer
comment et dans quelle mesure on peut rassembler ces données avec des courbes univer-
selles. On s’intéresse au comportement de la turbine Darrieus lorsque le confinement latéral
(ie la largeur du canal) varie. Les paramètres de calcul sont les mêmes que décrits dans la
partie 2. Le maillage du rotor est inchangé et celui du stator est reconstruit pour chaque
confinement. On définit le paramètre de confinement ǫ tel que :
W
ǫ= (3.15)
D
3.2. DESCRIPTION DU MODÈLE 109
ǫ prend cinq valeurs comprises entre 1.7 et 5.7. Le rapprochement des parois va avoir pour
effet une concentration du flux incident dans la turbine en empêchant l’expansion de la
veine fluide comme présenté pour l’écoulement libre. La figure 3.6 présente l’allure des
lignes de courant pour ǫ = 1.7 et 5.7. Elles montrent que la section amont de la veine fluide
est plus large dans le cas confiné.
On observe donc une augmentation forte du débit dans la turbine et par conséquent on
pressent une forte augmentation de la puissance pour une même vitesse amont. On re-
marque également une forte accélération du fluide dans la zone latérale. Cette accélération
a un effet sur la stabilité du sillage qui sera discuté dans le chapitre suivant. On reprend
les coefficients globaux définis au chapitre précédent par rapport à la vitesse amont. On
oubliera ici la notation m pour étudier une seule machine dans le canal :
ωR
λ= (3.16)
V∞
Cω
CP = (3.17)
ρRhV∞3
FX,Y
CX,Y = (3.18)
ρRhV∞2
La Figure 3.7 présente les coefficients de puissance pour différents confinements. On observe
une augmentation forte du CP optimal d’un facteur supérieur à 2, ainsi qu’un décalage
de la vitesse spécifique à l’optimum. En effet, à fort confinement, la vitesse vue par la
machine sera plus grande, elle doit donc tourner plus vite pour rester dans la zone de
transition. Un effet similaire est remarqué sur le CX ; les situations confinées montrent un
palier de saturation du coefficient de trainée. A partir de vitesses de rotation suffisamment
importantes, la vitesse débitante devient très faible et la force de trainée n’augmente plus.
On remarque que ce régime n’est pas atteint à λ = 4 pour la situation ǫ = 1.7.
Les courbes précédemment décrites permettent de connaître la performance de la machine
seule pour un confinement donné. Contrairement aux modèles analytiques comme celui de
Garett et Cumin [14] par exemple, le confinement ne doit pas intervenir comme paramètre
de notre modèle. L’augmentation de vitesse sera simulée par le calcul de l’écoulement.
Cela permet de prendre en compte également les effets de bord comme dans le cas d’une
configuration en ligne de quelques machines. L’idée décrite plus haut est donc de rapporter
la situation à des paramètres adimensionnels locaux, ie avec la vitesse interne en référence.
∗
Les Figures 3.8 a et b) montrent ces coefficients locaux CX et CP∗ en fonction du paramètre
d’avance local. Les courbes représentatives des différents confinements se rassemblent sur
∗
une courbe unique. CX montre une relation linéaire avec λ∗ avec une très bonne corrélation,
110 CHAPITRE 3. MÉTHODE D’ÉQUIVALENCE EMPIRIQUE
3.5
2.5
1.5
0.5
0
0.5 1 1.5 2 2.5 3 3.5 4 4.5
a)
1.2
0.8
0.6
0.4
0.2
−0.2
−0.4
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
b)
22
20
18
16
14
12
10
0
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
a)
−1
−2
−3
−4
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
b)
pour toutes les valeurs de ǫ. Le coefficient de puissance local CP∗ présente une courbe en
cloche. Les différentes configurations présentent des valeurs similaires de CP∗ pour les faibles
valeurs du paramètre d’avance local. Cet écart est dû au fait que l’effet de confinement n’est
pas équivalent à une augmentation de la vitesse de l’écoulement libre. Le forçage du courant
dans la veine a également un effet non linéaire sur la performance. Pour une même valeur
de λ∗ , les angles d’attaque d’une pale au cours de la rotation sont modifiés. Ainsi à grandes
valeurs de λ∗ on a un écart > 20% sur la valeur du CP∗ . On remarquera néanmoins que le
λ∗ associé à l’optimum du CP (global) se situe à des valeurs entre 4 et 5. Ce constat est
lié au résultat de l’analyse de Betz sur l’optimum de récupération de puissance, pour un
facteur d’induction constant. L’universalité d’une courbe pour CP∗ = f (λ∗ ) est donc assez
bonne dans la zone d’intérêt, c’est à dire les régions optimales et sub-optimales. On fait
donc le choix d’avoir une interpolation unique des ces courbes, avec comme seul paramètre
λ∗ .
Répartition spatiale Dans le cas des modèles de disque d’action pour des turbines à
flux axial, le terme source est en général imposé sur une ligne de cellule perpendiculaires
au flux. En raison de la complexité de l’écoulement dans le rotor de la Darrieus, il est peu
judicieux d’imposer la même répartition. Il est nécessaire d’obtenir un impact réaliste sur
l’écoulement, dans le sillage de la turbine, mais également une vitesse juste dans le rotor,
point clé du modèle. En étudiant l’évolution de la force de trainée imposée sur une pale au
cours de la rotation (Figure 3.9, on remarque deux régions de forte intensité). Une première
dans la zone amont θ = 70 − 120◦ , et la seconde à l’aval pour θ = 250 − 320◦ . L’amplitude
de la zone amont est environ deux fois supérieure à celle de la zone aval. Cette répartition
est globalement constante lorsque λ et ǫ varient. Il est donc judicieux d’appliquer une forme
similaire pour qX et qY . Une fonction d’espace arbitrairement choisie est imposée, que nous
appellerons la répartition I. Le terme source de quantité de mouvement est défini comme
le produit de deux fonctions d’espace :
1.6
1.4
1.2
0.8
0.6
0.4
0.2
−0.2
−50 0 50 100 150 200 250 300 350 400
Figure 3.9 – Coefficient de trainée sur une pale de la turbine Darrieus en fonction de la
position angulaire au point de fonctionnement optimal pour un blocage de ǫ = 4 ; l’origine
de θ est prise lorsque la pale remonte le courant
rotation pour que la valeur de λ calculée avec la vitesse en entrée du domaine soit égale à
2. On remarque que la vitesse apparaît plus faible dans la zone du disque dans le cas de la
répartition réaliste. En moyenne sur le disque, la vitesse est 30% plus importante. On ob-
tient également une surestimation du coefficient de puissance avec une valeur de 0,69 contre
0,49. Les figures a) et b) sont issues des mêmes calculs et représentent le champ d’énergie
cinétique turbulente. On remarque une zone de forte production turbulente convectée à
l’aval créée par la répartition complexe. Le gradient de force volumique engendre un gra-
dient de vitesse et ainsi une zone de fort cisaillement. La production de k participe ensuite
au phénomène de mélange dans le sillage, où l’on observe une réalimentation plus rapide
avec notre répartition que dans le cas du disque uniforme. L’influence de la répartition de
source de quantité de mouvement et d’énergie cinétique turbulente est discutée au chapitre
suivant, afin de comprendre l’effet sur le sillage reproduit.
a) b)
Figure 3.10 – Champ de vitesse autour de la turbine, obtenus avec le modèle pour la
répartition I (a) et pour une répartition uniforme (b)
et une condition de sortie standard de Code_Saturne (ie une condition de Neumann ho-
∂P
mogène pour toutes les grandeurs transportées et ∂n∂τ i
= 0, i ∈ [1, 2], avec n le vecteur
normal à la face de sortie et τi deux vecteurs orthogonaux entre eux dans le plan de la face
de sortie.
Le choix du modèle de turbulence a un impact sur la stabilité de l’écoulement cisaillé. Il
a donc également un impact sur la convergence du schéma itératif pour le calcul du point
de fonctionnement de la turbine. On remarque que quelque soit le schéma en temps utilisé
(stationnaire ou instationnaire), le champ de vitesse est instable dans le sillage de la turbine
ce qui a un impact sur la vitesse locale et donc sur le CP calculé. Le nombre de mailles
nécessaires pour une bonne reproduction de la répartition des forces est également étudié.
Le tableau 3.1 récapitule les différents tests effectués avec une turbine dans le domaine de
la veine d’essai à λ =2.
Table 3.1 – Test de convergence de la méthode pour différents maillages, schémas tempo-
rel, ou modèles de turbulence
1.2
maillage1
1.1 maillage2
maillage3
1 maillage4
maillage5
0.9
0.8
0.7
CP
0.6
0.5
0.4
0.3
0.2
0.1
0
0 100 200 300 400 500 600 700 800 900 1 000
iteration
taille de maille inférieure à R/8 qui correspond au maillage 4 dans la zone de l’hydrolienne
pour la suite des calculs. Les variations et instabilités sont expliquées au chapitre suivant
par la faible production en modèle v 2 − f et k − ǫ-PL par rapport au kω − SST ou au
k − ǫ. Les forts cisaillements créés par les forces volumiques entraînent une instabilité de
l’écoulement, voir des recirculations dont l’effet est très mauvais pour la convergence du
schéma itératif.
0.9
0.8
0.7
0.6
CP
0.5
0.4
0.3
0.2
0.1
0
0 100 200 300 400 500 600 700 800 900 1 000
iteration
(a) (b)
3.3 Validation
Le travail effectué sur la méthode de maillage glissant de Code_Saturne et les bons ré-
sultats obtenus apportent des résultats fiables sur les performances de la turbine dans des
configurations variées. Une partie du développement a également permis la simulation de
plusieurs machines dans un même domaine. Les données issues de ces calculs servent à la
validation de la méthode décrite dans ce chapitre. Les résultats des deux simulations sont
comparés dans le cas d’une machine simple puis pour des parcs de quelques machines. La
fiabilité de la méthode nous permettra de l’utiliser pour la compréhension des phénomènes
généraux dans les parcs.
Le maillage utilisé pour le calcul avec le modèle simplifié est conforme et orthogonal avec
une taille de maille uniforme de 1cm2 . La turbine est ainsi représentée sur environ 300
mailles. Les gradients de vitesse créés dans le sillage de la perte de charge engendrent
un production turbulente, et des phénomènes d’instabilité conduisent à l’utilisation d’un
calcul instationnaire. Le pas de temps est de l’ordre de la milliseconde pour un nombre de
Courant de l’ordre de l’unité. Le calcul est mené sur 1000 itérations, et les résultats sont
moyennés sur les 500 dernières.
Les résultats pour les 3 valeurs du confinement sont présentés sur la figure 3.13. Les courbes
de référence sont en pointillé et les résultats du modèle sont représentés par des points.
Les coefficients de puissance donnés par le modèle suivent fidèlement les valeurs attendues
pour ce cas simple. De façon plus exacte, on remarque que la concordance est la meilleure
pour les confinements faibles (ǫ = 5.7 et 4). Pour un blocage plus fort, un décalage des
courbes est visible. Le modèle prédit l’optimum du CP pour des λ∗ légèrement plus élevés
que ce que l’on attend.
Dans ce chapitre, nous concentrons notre attention sur la prédiction des performances.
Ainsi nous ne prendrons que la valeur du CP comme critère de validation. La restitution
du sillage est examinée au chapitre suivant.
118 CHAPITRE 3. MÉTHODE D’ÉQUIVALENCE EMPIRIQUE
1.2
0.8
0.6
0.4
0.2
−0.2
−0.4
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
Machine 1
Wall
1D
Stream−wise Velocity
2D
1.5
U/Uref
Machine 2
11.5 D
0.5
Inlet Outlet
equivalence model
0
full CFD
Transverse Velocity
0.2
−0.2
8D
V
−0.4
dèle. La partie supérieure qui voit la présence des turbines dispose d’un maillage structuré.
Les mailles sont du même ordre de grandeur que pour le cas test précédent. Le reste du
domaine est maillé plus grossièrement, afin d’économiser en temps de calcul.
Le tableau 3.2 présente les résultats des simulations avec les deux méthodes. La dernière
colonne donne l’erreur relative sur le coefficient de puissance moyen des deux machines.
On remarque un bonne concordance des résultats pour les valeurs du paramètre d’avance
jusqu’à λ = 2, 5. L’optimum de puissance pour les deux machines semble être situé aux
alentours de cette valeur, ce qui s’explique par l’effet de confinement créé par la présente
disposition. Les erreurs plus importantes sont remarquées pour les points de fonctionnement
au-delà de l’optimum (λ = 3), pour lesquels le modèle par perte de charge surestime la
valeur de CP .
120 CHAPITRE 3. MÉTHODE D’ÉQUIVALENCE EMPIRIQUE
Table 3.2 – Résumé des résultats du cas de validation asymétrique à deux turbines ; la
dernière colonne rapporte l’erreur relative l’erreur sur la valeur moyenne du coefficient de
puissance
Points clés
• Ce chapitre présente la méthodologie de construction d’un modèle simplifié de turbine
à partir d’une paramétrisation judicieuse des données empiriques (issues d’un calcul
complet).
• Le modèle décrit se place dans une même classe d’utilisation que les modèles
RANS+BE [26, 19], ie il peut servir pour l’optimisation de placement au sein de
grappe de machines dans un parc. L’utilisation de données empiriques sur une plage
étendue de points de fonctionnement permet un calcul moins lourd de l’action de la
turbine dans la simulation d’écoulement.
• Le fonctionnement du modèle est validé par comparaison avec les résultats de calcul
en géométrie complète pour une turbine seule dans un canal.
• La performance est calculée correctement pour un confinement faible. On observe
une surestimation au delà de l’optimum pour le confinement plus fort (à partir de
ǫ = 2, 7).
• La comparaison du sillage obtenu par l’expérience et les différentes simulation
montrent une réalimentation trop rapide avec le modèle simple. Il semble que la
production turbulente surestimée produise un mélange important qui contribue à la
diffusion de quantité de mouvement.
• La répartition des forces dans la machine a un effet important sur la production
d’énergie cinétique turbulente et par conséquent la stabilité du schéma itératif.
BIBLIOGRAPHIE 121
Bibliographie
[1] A.-M. Georgescu, S.-C. Georgescu, C. I. Cosoiu, and N. Alboiu, “Efficiency of marine
hydropower farms consisting of multiple vertical axis cross-flow turbines,” Internatio-
nal Journal of Fluid Machinery and Systems, 2011.
[3] T. Stallard, R. Collings, T. Feng, and J. Whelan, “Interactions between tidal turbine
wakes : experimental study of a group of three-bladed rotors,” Philosophical Tran-
sactions of the Royal Society A : Mathematical, Physical and Engineering Sciences,
vol. 371, no. 1985, p. 20120159, 2013.
[4] A. Bahaj and L. Myers, “Shaping array design of marine current energy converters
through scaled experimental analysis,” Energy, vol. 59, pp. 83–94, 2013.
[5] K. Bergey, “The lanchester-betz limit (energy conversion efficiency factor for wind-
mills),” Journal of Energy, vol. 3, no. 6, pp. 382–384, 1979.
[6] F. W. Lanchester, “A contribution to the theory of propulsion and the screw propeller,”
Journal of the American Society for Naval Engineers, vol. 27, no. 2, pp. 509–510, 1915.
[9] C. Garrett and P. Cummins, “The efficiency of a turbine in a tidal channel,” Journal
of fluid mechanics, vol. 588, pp. 243–251, 2007.
[10] J. Whelan, J. Graham, and J. Peiro, “A free-surface and blockage correction for tidal
turbines,” Journal of Fluid Mechanics, vol. 624, pp. 281–291, 2009.
[11] T. Nishino and R. H. Willden, “The efficiency of an array of tidal turbines partially
blocking a wide channel,” Journal of Fluid Mechanics, vol. 708, pp. 596–606, 2012.
[12] T. Nishino and R. H. Willden, “Two-scale dynamics of flow past a partial cross-stream
array of tidal turbines,” Journal of Fluid Mechanics, vol. 730, pp. 220–244, 2013.
122 CHAPITRE 3. MÉTHODE D’ÉQUIVALENCE EMPIRIQUE
[13] R. Vennell, “Exceeding the betz limit with tidal turbines,” Renewable Energy, vol. 55,
pp. 277–285, 2013.
[14] C. Garrett and P. Cummins, “The power potential of tidal currents in channels,”
Proceedings of the Royal Society A : Mathematical, Physical and Engineering Science,
vol. 461, no. 2060, pp. 2563–2572, 2005.
[15] R. Vennell, “Estimating the power potential of tidal currents and the impact of power
extraction on flow speeds,” Renewable Energy, vol. 36, no. 12, pp. 3558–3565, 2011.
[16] M. Islam, D. S.-K. Ting, and A. Fartaj, “Aerodynamic models for darrieus-type
straight-bladed vertical axis wind turbines,” Renewable and Sustainable Energy Re-
views, vol. 12, no. 4, pp. 1087–1109, 2008.
[17] G. I. Gretton, Hydrodynamic analysis of a vertical axis tidal current turbine. PhD
thesis, 2009.
[19] S. Antheaume, T. Maître, and J.-L. Achard, “Hydraulic darrieus turbines efficiency for
free fluid flow conditions versus power farms conditions,” Renewable Energy, vol. 33,
no. 10, pp. 2186–2198, 2008.
[20] M. Harrison, W. Batten, L. Myers, and A. Bahaj, “Comparison between cfd simulations
and experiments for predicting the far wake of horizontal axis tidal turbines,” IET
renewable power generation, vol. 4, no. 6, pp. 613–627, 2010.
[22] R. Malki, A. Williams, T. Croft, M. Togneri, and I. Masters, “A coupled blade element
momentum–computational fluid dynamics model for evaluating tidal stream turbine
performance,” Applied Mathematical Modelling, vol. 37, no. 5, pp. 3006–3020, 2013.
[23] G. Bai, J. Li, P. Fan, and G. Li, “Numerical investigations of the effects of different ar-
rays on power extractions of horizontal axis tidal current turbines,” Renewable Energy,
vol. 53, pp. 180–186, 2013.
BIBLIOGRAPHIE 123
[24] A. Bahaj, A. Molland, J. Chaplin, and W. Batten, “Power and thrust measurements
of marine current turbines under various hydrodynamic flow conditions in a cavitation
tunnel and a towing tank,” Renewable energy, vol. 32, no. 3, pp. 407–426, 2007.
[25] T. Roc, D. C. Conley, and D. Greaves, “Methodology for tidal turbine representation
in ocean circulation model,” Renewable Energy, vol. 51, pp. 448–464, 2013.
[26] W. M. Batten, M. Harrison, and A. Bahaj, “Accuracy of the actuator disc-rans ap-
proach for predicting the performance and wake of tidal turbines,” Philosophical Tran-
sactions of the Royal Society A : Mathematical, Physical and Engineering Sciences,
vol. 371, no. 1985, 2013.
[27] L. Myers and A. Bahaj, “Near wake properties of horizontal axis marine current tur-
bines,” in Proceedings of the 8th European Wave and Tidal Energy Conference, pp. 558–
565, Uppsala, Sweden, 2009.
[28] V. Aumelas, Modélisation des hydroliennes à axe vertical libres ou carénées : développe-
ment d’un moyen expérimental et d’un moyen numérique pour l’étude de la cavitation.
PhD thesis, Grenoble-INP, 2011.
124 CHAPITRE 3. MÉTHODE D’ÉQUIVALENCE EMPIRIQUE
Chapitre 4
Le sillage des turbines hydrocinétiques, qu’elles fonctionnent dans l’air ou bien dans l’eau,
est caractérisé par deux phénomènes principaux : la région à l’aval de la machine voit
un déficit de quantité de mouvement par rapport à l’écoulement libre, le second concerne
l’augmentation du niveau de turbulence [1]. Ces deux effets ont un impact important sur
la puissance récupérable dans les parcs et sur la sollicitation mécanique des structures. La
compréhension des sillages est donc un point crucial pour la détermination des dispositions
de machines dans un parc. Ces interactions entre machines nécessitent un soin particulier,
pour comprendre l’effet des différents paramètres sur la performance. A titre d’exemple,
une augmentation de l’intensité turbulente a des effets mal maîtrisés sur l’efficacité de la
turbine, mais favorise la récupération de la vitesse dans le sillage [2]. Pour autant, les
fluctuations de grandes échelles créent des chargements variables et une importante fatigue
mécanique. On ne peut donc pas faire l’économie d’une connaissance globale et structurelle
de l’écoulement qui attaquera une éventuelle deuxième rangée de turbines.
Cette problématique du sillage fait l’objet de nombreuses études, à la fois analytiques,
expérimentales et numériques. Les travaux de synthèse bibliographique de Crespo, Veermer
et Sanderse [1, 3, 4] sont des sources d’information de grande valeur. Ils sont mis en évidence
plusieurs phénomènes physiques déterminants pour l’établissement, puis la dissipation du
sillage.
Plusieurs phénomènes physique principaux rentrent en compte pour l’établissement puis la
dissipation du sillage. Les principaux phénomènes se résument comme suit :
125
126 CHAPITRE 4. ETUDE DU SILLAGE DE LA TURBINE ACHARD
Dans ce chapitre, une revue des connaissances concernant le sillage des éoliennes et hydro-
liennes est exposée, ainsi qu’une analyse des méthodes de simulation employées. Elle met
en avant le peu de connaissance sur le sillage des turbines à flux transverse. Nous tentons
de comprendre les différences entre les deux types de turbines pour savoir si les modèles
proposés pour les turbines classiques sont utilisables. Des mesures expérimentales de vi-
tesses dans le sillage d’un modèle réduit d’hydrolienne dans la boucle hydrodynamique du
LEGI permettent d’identifier les phénomènes dominants dans le sillage proche de la tur-
bine. Enfin, les modèles RANS en géométrie complète et par le modèle d’équivalence sont
examinés pour leur réponse en fonction du modèle de turbulence utilisé et de paramètres
d’entrée. L’objectif général de cette section est de déterminer les capacités du modèle sim-
plifié à reproduire le sillage de la turbine afin de pouvoir simuler des parcs de deuxième
génération. Le déficit de vitesse et l’augmentation de la turbulence à différentes échelles
doivent être maîtrisés, et la pertinence des modèles RANS pour cette région sera discutée.
Figure 4.1 – Visualisation des tourbillons de bout de pale d’une éolienne par un panache
de fumée. A gauche, la fumée est émise à l’amont du rotor, à droite un dispositif émet la
fumée au bout de la pale (tirés respectivement de Hand [6] et Alfredsson [7])
Les éoliennes à flux axial Le sillage de l’éolienne, regardé de façon globale, est large-
ment dépendant de la force de trainée appliquée par le fluide sur la turbine. L’extraction
d’énergie par la turbine, provoque un ralentissement dans la zone du rotor et un contour-
nement de l’écoulement. Le déficit de vitesse est donc relié à la trainée, qui peut être
vue comme la perméabilité du rotor. Plus il est bloquant (typiquement plus le paramètre
d’avance est fort), plus le déficit de vitesse à l’aval immédiat de la turbine sera important.
La figure 4.2 qui présente un profil de vitesse adimensionnelle dans le sillage à X/D=1,67
pour plusieurs paramètres d’avance [8], illustre ce phénomène.
Le sillage d’une turbine est habituellement séparé en deux zones principales. La zone du
sillage proche, située juste à l’aval de la turbine, est le siège du plus fort déficit en vitesse.
Dans cette région le sillage est dépendant de la géométrie de la turbine. Les tourbillons créés
par les pales y sont convectés et la vitesse de réalimentation est faible. Des visualisations
des tourbillons de bout de pales par des panaches de fumée sont visibles sur la figure 4.1.
Dans les deux cas des tourbillons sont formés dans la direction ortho-radiale par l’extrémité
de la pale et sont transportés par l’écoulement.
La zone du sillage lointain en aval du sillage proche s’étend jusqu’à la réalimentation de
la vitesse axiale. Le mélange est affecté par deux phénomènes qui sont la convection et
la diffusion turbulente. Le sillage s’étend latéralement jusqu’à retrouver la vitesse initiale.
Dans le cas des turbines axiales, les tourbillons lâchés par les extrémités de pale sont
convectés vers l’aval (figure 4.1). Ils forment deux couches fines de cisaillement, qui sont la
frontière entre le sillage et l’écoulement dévié par la turbine. La diffusion turbulente tend
à élargir la couche de cisaillement en s’éloignant du rotor. La frontière entre sillage proche
128 CHAPITRE 4. ETUDE DU SILLAGE DE LA TURBINE ACHARD
=0
U/U0
=5,0
=6,6
=8,5
r/R
a) b)
Figure 4.2 – Gauche, déficit de vitesse d’un éolienne à flux axial, expérience de Vermeulen
[8], Droite, vue schématique de l’écoulement dans le sillage (d’après Bahaj et Myers [9])
et lointain est définie comme la distance à partir de laquelle ces couches de cisaillement
atteignent l’axe du sillage (voir figure 4.2.b).
De nombreuses études expérimentales permettent de comprendre finement la structure du
sillage et l’influence relative des phénomènes physiques sur le mélange. Zhang et al. ana-
lysent le sillage proche d’une éolienne dans un écoulement de couche limite atmosphérique
turbulente [10]. L’intensité turbulente dans le tunnel est d’environ 2%. Le fluide dans le
sillage est étudié par mesures PIV dans des plans transverses à X/D = 1, 2, 3, 5 à l’aval du
rotor et dans le plan vertical à l’axe de la machine. Ils montrent que les tourbillons de bout
de pale ont une forte empreinte sur le sillage proche jusqu’à une distance de 2-3 diamètres
à l’aval. Ces structures sont mises en évidence par la présence d’un pic de cohérence tempo-
relle dans l’analyse spectrale de la vitesse longitudinale. En raison de l’asymétrie verticale
de la couche limite, les tourbillons sont plus marqués et l’intensité turbulente plus forte
dans la partie haute du sillage. Ils ont mis en évidence un élargissement du sillage (distance
entre les zones cisaillées) jusqu’à X/D = 1 et montrent que les intensités turbulentes sur
les 3 composantes de l’espace ont des valeurs comparables dans le sillage proche. A partir
4.1. CONSIDÉRATIONS SUR LE SILLAGE DES HYDROLIENNES 129
Kang et al. [14] proposent une analyse numérique des mesures de Chamorro en modèle de
turbulence LES, dans laquelle la géométrie complète de la turbine est décrite à l’aide de
la méthode des frontières immergées. Ils obtiennent une très bonne concordance lors de la
confrontation aux résultats expérimentaux. La même importance des structures cohérentes
est retrouvée. Ils estiment que l’élargissement du sillage et son interaction avec l’écoulement
latéral provoquent la destruction de ces structures après 4 diamètres à l’aval du rotor. Cette
distance correspond à la zone de déstabilisation à grande échelle du sillage caractérisée par
une forte hausse de l’énergie cinétique turbulente (TKE). Ce phénomène d’instabilité à
grande échelle du sillage est appelé meandering. Deux explications principales sont données :
1. Le sillage fonctionnerait comme un amplificateur des fluctuations de la vitesse amont.
Il provient des variations dans la direction du vent. D’après Larsen [13], le sillage
agit comme un traceur passif des grandes structures atmosphériques, illustré sur la
figure 4.3.
2. Ce phénomène est similaire aux instabilités de grandes tailles visibles derrière un
corps bloquant dans un écoulement. Medici [12] étudie le meandering dans le sillage
130 CHAPITRE 4. ETUDE DU SILLAGE DE LA TURBINE ACHARD
d’une éolienne en utilisant une soufflerie à la turbulence contrôlée (figure 4.4.a). Cette
figure représente la corrélation temporelle entre deux points de mesures de la vitesse
(uX à gauche, uY à droite) situés dans le sillage. On voit deux fréquences de struc-
tures cohérentes, dus respectivement aux lâchers des extrémités de pales, et à des
fluctuations grande échelle. Le phénomène intervient dans les situations avec et sans
augmentation de la turbulence. Ils constatent de plus qu’il apparaît à partir d’une
certaine valeur minimum du coefficient de traînée qui augmente avec la vitesse de ro-
tation relative [12],(figure 4.4.b). L’instabilité du sillage est retrouvée avec des disques
poreux jusqu’à 42% de porosité [15].
a) b)
Figure 4.4 – a) Corrélation temporelle entre la vitesse mesurée en deux points situés à
X/D=1 et à distance Y variable dans le sillage d’une éolienne. Les flèches noires indiquent
le passage d’une pale, et les flèches blanches les fluctuations à grande échelle. b) Coefficient
de trainée pour différents paramètres de l’éolienne. Les points en gras indiquent la présence
du battement du sillage. Extrait de Medici et al. [16, 12]
de mesure est menée pour des intensités turbulentes amont de It = 8% et It = 25%. Une
récupération de la vitesse est atteinte après 4D seulement pour le cas hautement turbulent
contre plus de 10D pour le cas moyen.
Mycek et al. étendent la même étude en s’intéressant tout d’abord à l’effet de la turbulence
sur la performance [17]. La puissance moyenne récupérée par le modèle réduit de turbine
sous It = 3 ou 15% est peu différente quoique légèrement inférieure à forte intensité turbu-
lente. On peut noter toutefois que les fluctuations du CP sont fortement augmentées, ce qui
sera préjudiciable d’un point de vue mécanique et électrique. La situation est également
envisagée en positionnant une deuxième turbine à une distance variable derrière la première
[18, 19]. Ils établissent ainsi qu’il est avantageux de fonctionner en régime turbulent pour
le positionnement d’hydrolienne en parc afin de réduire la distance entre les lignes. Pour
une distance a/D = 6 et It = 3% ils constatent une chute de CP de 75% contre seulement
20% à It = 15%.
Dans la même optique, Stelzenmuller et Aliseda étudient trois modèles réduits dans un
canal [20]. Grâce à la production turbulente de la première turbine, la seconde fonctionne
à un niveau de turbulence ambiant élevé, ce qui accélère la récupération de la vitesse axiale
entre la turbine 2 et 3. De façon surprenante, une puissance similaire est obtenue pour les
turbines 1 et 3 à certains espacements.
La fiabilité des modèles simplifiés dans le sillage : RANS Les études numériques
sur des écoulements atmosphériques de canopées utilisent des techniques de termes de
quantité de mouvement pour représenter l’interaction du vent sur la végétation, comme
des forces de trainée. Sanz [27] propose un modification du modèle k − ǫ par l’adjonction
de termes sources dans les équations de conservation de k et de ǫ afin de prendre en compte
deux effets :
2. le transfert d’énergie turbulente des grandes échelles vers les petites appelé également
short-circuiting of tubulent cascad [28].
1
SU = − C x U 2 (4.2)
2
1
Sk = Cx (βp U 3 − βd U k ) (4.3)
2 | {z } | {z }
dissip non−linear
1 ǫ
Sǫ = Cx (Cǫ5 βp U 3 − Cǫ5 βd U ǫ) (4.4)
2 k
SU est donnée par la définition classique du coefficient de trainée, proportionnel au carré
de la vitesse. En ce qui concerne Sk , le premier terme, avec βp ∈ [0, 1], est relatif à l’énergie
cinétique perdue par le fluide transformée en énergie turbulente. Le deuxième correspond
au transfert d’énergie entre échelles qui est modélisée comme une diminution de l’énergie
cinétique turbulente. Réthoré et al. proposent d’adapter le modèle de Sanz aux éoliennes
[29]. Cette adaptation passe par une estimation de βp par la théorie du disque d’action.
Cette constante est reliée à l’énergie cinétique perdue par le fluide sous forme de turbulence,
et prise égale à 1 par Sanz. En réalité, ce terme doit également prendre en compte les pertes
sous forme de chaleur par dissipation visqueuse et l’énergie récupérée mécaniquement par
le générateur, donnant : P = Pelec + Pturb + Pchaleur . Ils estiment ainsi le premier membre
de Sk :
1 Cchaleur + CP 3
P = ρACx (1 − )U (4.5)
2 4a(1 − a)2
le terme 4a(1 − a) est le coefficient de trainée à l’optimum dans la théorie de Betz. En
prenant CP = 0, 4 et CX = 0, 8, ils trouvent βp = 0, 05. βd = 2/3 est donné analytiquement
en appliquant la décomposition des Reynolds sur les termes sources. Les auteurs trouvent
4.2. MESURES EN SILLAGE PROCHE DE LA DARRIEUS 135
une amélioration de la longueur de réalimentation donnée par leur modèle et des résultats
plus proches du calcul LES.
El Kasmi propose un modèle différent pour tenir compte de l’interaction non linéaire
d’échelle [30] en modèle k − ǫ. Il utilise le modèle de disque d’action uniforme pour si-
muler une éolienne. Deux modifications sont effectuées dans l’équation de conservation du
taux de déformation :
1. Une nouvelle échelle de temps est ajoutée pour améliorer la réponse du taux de
dissipation au taux de déformation moyen, qui ajoute un terme de production dans
l’équation de ǫ.
2. L’ajout d’un terme source Pǫ , qui est appliqué uniquement dans une zone du rotor
(eq :(4.6)). La zone s’étend pour leur calcul à 2,5D autour de la machine. D’après
l’auteur, le choix de cette zone a une forte influence sur les résultats mais il n’explique
pas laquelle.
P2
Pǫ = Cǫ4 t (4.6)
ρk
Les résultats du modèle modifié sont comparés avec les mesures de Taylor, 1985, pour
une intensité turbulente de 11%, et différentes chargements, CX = 0, 67; 0, 77; 0, 82. La
restitution du déficit de vitesse est nettement améliorée par rapport au k − ǫ classique de
X/D=2,5 à 7,5 pour CX = 0, 67 mais le modèle surestime le déficit dans le sillage lointain
avec des coefficients de trainée plus importants. L’auteur ne compare pas l’influence relative
des deux modifications.
Principe On utilise un montage de LDV à faisceaux croisés. Le point de mesure est éclairé
par deux faisceaux cohérents issus du même laser. Le croisement de deux faisceaux crée des
franges d’interférence, dans une région restreinte, que l’on appellera volume de mesure. Elles
sont dues à la différence de marche dans le chemin optique de chacun des faisceaux, de leur
émission jusqu’au point de mesure. Une particule traversant cette région voit une alternance
de zones éclairées et sombre, espacées régulièrement. Elle émet des flashs lumineux à une
fréquence proportionnelle à sa vitesse de déplacement dans la direction orthogonale aux
franges. Le principe général de fonctionnement est illustré sur la figure 4.5. La distance
interfrange dépend de l’angle entre les faisceaux 2θ et la longueur d’onde du faisceaux :
λ
df = (4.7)
2sinθ
Ce dispositif présente plusieurs avantages :
• La distance interfrange (df ) est indépendante du milieu de diffusion. Dans notre cas
le laser émet dans l’air et la mesure de vitesse est effectuée dans l’eau. Par la loi de
λair
Descartes nair λair = neau λeau et nair sinθair = neau sinθeau et ainsi : df = 21 sin(θ air )
=
1 λeau
2 sin(θeau )
.
Ensemencement Les particules choisies pour les mesures dans la veine d’essais du LEGI
sont des billes de verre creuses de 10µm recouvertes d’une couche d’argent réfléchissante.
La quantité de particules est difficilement maîtrisable du fait de l’important volume d’eau,
et du phénomène de sédimentation qui intervient dans les zones de stabilisation. Il s’agit de
trouver un compromis qui donne le meilleur taux d’acquisition sans avoir de superposition
de particules dans le volume de mesure. La limitation est donc imposée par le temps de
transit de la particule dans le volume. Dans notre cas on a un volume de mesure de 2.10−4 m
et une vitesse moyenne de l’écoulement de 2, 3m/s. La vitesse d’acquisition maximale est
donc fmax = 10kHz. Suivant la zone de l’écoulement mesurée, on atteint effectivement des
fréquences de l’ordre de 2-3 kHz.
Figure 4.6 – Le dispositif de mesure LDV est placée sous la veine d’essai, supporté par
un chariot de déplacement automatique
4.2. MESURES EN SILLAGE PROCHE DE LA DARRIEUS 139
0.6 0.2
0.5
0.4 0.15
0.3
0.2 0.1
0.1
0 0.05
−0.1
−0.2 0
−0.3
−0.4 −0.05
−0.5
−0.6 −0.1
−0.7
−0.8
T T T −0.15
11 11.05 11.1 11.15 11.2 11.25 11.3 11.35 11.4 0 50 100 150 200 250 300 350 400
Temps aquisition
a ) v (m/s) b ) Ṽ (m/s)
A partir de cette décomposition, on calcule les moyennes des termes diagonaux du tenseur
de Reynolds, u′2 et v ′2 . L’énergie cinétique turbulente k est donnée en 2D par la relation
suivante :
1
k = (u′2 + v ′2 ) (4.9)
2
Afin de calculer le moment croisé u′ v ′ (terme non diagonal du tenseur de Reynolds), les évé-
nements enregistrés en non-coïncidence doivent être reconstruits a posteriori. Pour chaque
événement sur U, on identifie la correspondance sur V. On fixe pour cela un critère sur
la date d’arrivée du signal lumineux de ±3.10−5 s ce qui est inférieur au temps de transit
minimum d’une particule dans le volume de contrôle. La reconstruction permet de garder
en moyenne la moitié des événements.
La résolution spatiale de la grille de mesure permet de calculer le taux de rotation et de
déformation du champ moyen, définis de la façon suivante :
∂U ∂V
S12 = ( + ) (4.10)
∂Y ∂X
4.2. MESURES EN SILLAGE PROCHE DE LA DARRIEUS 141
∂V ∂U
Ω12 = ( − ) (4.11)
∂X ∂Y
Le taux de production d’énergie cinétique turbulente k est défini comme le produit du
tenseur de Reynolds, et du gradient du champ moyen, par l’équation de conservation de
l’énergie cinétique turbulente (voir chapitre 2). C’est donc la contribution de quatre termes.
Pour un écoulement quasi-stationnaire comme dans notre cas, la contribution principale
provient de u′ v ′ ∂U
∂y
.
∂Ui
P rod = u′i u′j (4.12)
∂xj
Une seconde décomposition est appliquée. Elle est plus proche de la signification de la
décomposition de Reynolds et permet la comparaison avec les résultats du calcul RANS.
Elle est utile pour évaluer le poids des fluctuations périodiques sur l’énergie turbulente
totale. On définit alors la vitesse de perturbation Ũ .
u(t) = U + Ũ (t) + u′′ (t) (4.13)
Pour obtenir Ũ , on effectue une moyenne de phase du signal temporel en chaque point,
échantillonnée sur 30 intervalles. L’acquisition n’est pas synchronisée sur la position de la
turbine, en raison de l’échantillonnage aléatoire du signal. On décompose le signal brut en N
tronçons d’une durée égale à la période de rotation T (figure 4.7). Chaque tronçon est divisé
en P intervalles (ici P = 30). La valeur Ũ (i) est calculée en prenant la moyenne de tous
les événements appartenant au j e me intervalle de chaque tronçon. Cette décomposition
a l’avantage d’être simple et rapide pour un signal dont la périodicité est connue, par
rapport à un filtrage fréquentiel. On regarde seulement la fluctuation déterministe associée
au passage des pales et on élimine les autres. On peut voir sur le spectre de la vitesse
calculé dans la zone cisaillée pour deux points à différentes distances à l’aval (figure 4.8),
que la fréquence associée aux tourbillons évolue.
On peut alors calculer la variance de ce signal qui est représentative de l’énergie transportée
par les grandes structures vorticitaires, résolues en RANS. On note K̃ la moyenne de la
variance sur chaque composante. Pour avoir accès à l’énergie cinétique turbulente seule, on
utilise la formule (4.14). La contribution du terme de corrélation des deux fluctuations est
inférieure de 2 ordres de grandeur par rapport à k.
sens anti-horaire, celui pris dans les simulations). Le champ de vitesse (4.9.a) montre un
début de réalimentation rapide (à partir de X/D=1) marqué par l’élargissement des zones
cisaillées (4.9.b). Le déficit de vitesse le plus fort est obtenu aux environ de X/D = 1.
Deux zones sont marquées par une forte énergie cinétique turbulente (4.9.c et d) . Elles
correspondent aux couches de cisaillement, mais également à la trajectoire des tourbillons
lâchés périodiquement.
Le sillage global est fortement asymétrique. En bas, l’intensité de turbulence marque un
pic aux environs de X/D = 1, 5, puis diminue. En haut, elle démarre dans cette zone.
La figure 4.9.f présente l’énergie de structures cohérentes. Elle permet de comprendre la
dynamique des deux parties. On voit que dans la partie basse, des structures de forte
amplitude sont lâchées, et que l’énergie associée décroit rapidement. En haut, le passage de
structures intervient également, mais avec une intensité moindre. Ce phénomène est plus
fortement marqué sur la vitesse transverse v que sur la vitesse axiale u.
Le déséquilibre est attribuable à la différence d’intensité entre les deux lâchers de tourbillons
dans la turbine. La dynamique de production dans le sillage proche semble donc être
pilotée par la dislocation des structures cohérentes jusqu’à X/D=1,5-2 qui transmettent
leur énergie vers les plus petites échelles d’agitation. Ensuite le cisaillement dû au déficit
de vitesse devient prépondérant et on retrouve la symétrie au delà de X/D=2, avec un
niveau de production turbulente équivalent. Dans l’axe de la turbine le niveau d’intensité
turbulente est faible. L’augmentation de l’agitation turbulente à l’aval de la machine est
donc due en partie à la création des tourbillons, mais surtout au cisaillement entre le sillage
et l’extérieur (figure 4.9 c et d).
L’aspect instationnaire du sillage proche disparaît rapidement et joue seulement un rôle
4.2. MESURES EN SILLAGE PROCHE DE LA DARRIEUS 143
dans l’asymétrie du sillage au démarrage. On remarque que la force globale sur la machine
dans la direction Y entraîne un décalage de tout le sillage vers le bas. On peut conjecturer
que cette déviation peut amener des instabilités à plus grandes échelles.
Figure 4.9 – Carte des grandeurs cinétiques dans le plan horizontal XY au quart des
pales.
144 CHAPITRE 4. ETUDE DU SILLAGE DE LA TURBINE ACHARD
Les conditions aux limites Dans un premier temps on veut s’affranchir de la produc-
tion turbulente induite par la présence de la turbine. On simule pour ce faire un écoulement
cisaillé dans un domaine de 30D*4D. On peut ainsi évaluer le comportement de différents
modèles RANS en présence d’un écoulement cisaillé. La condition d’entrée est un profil
défini par la condition (4.15) :
πy
U (y) = 0, 5 + 2sin2 ( ) (4.15)
4D
Pour les modèles RANS dans Code_Saturne, le taux de turbulence dans le fluide entrant
est réglé en imposant une valeur de l’intensité turbulente It et une valeur du diamètre hy-
draulique DH . Les valeurs de k et ǫ imposées à l’entrée sont déterminées selon les formules,
avec Cµ = 0, 09.
2
k = (It U0 )2 (4.16)
3
3/4
Cµ k 2/3
ǫ= (4.17)
0, 07DH
Le taux de dissipation est quant à lui inversement proportionnel à la longueur de mélange.
Plus la longueur de mélange est prise petite, plus la dissipation de l’énergie turbulente va
être importante. La distance entre l’entrée du domaine et la/les régions d’intérêt a donc
une influence sur la simulation, et notamment pour des valeurs de DH petite.
Ce profil donne des pentes de la vitesse axiale U dans la direction transverse moins impor-
tante que ce qui est constaté avec une turbine de façon expérimentale ou numérique. De ce
fait, on ne constate aucun phénomène instable. Néanmoins, il permet de regarder le phéno-
4.4. ETUDE NUMÉRIQUE DU SILLAGE EN MAILLAGE ROTATIF 145
1 5e−01
4.5e−01
0.9
4e−01
0.8 3.5e−01
3e−01
0.7
2.5e−01
0.6
2e−01
0.5 1.5e−01
1e−01
0.4
5e−02
0.3 0e00
0 5 10 15 20 25 30 35 0 5 10 15 20 25 30 35
(a) (b)
Figure 4.10 – Vitesse minimale dans l’axe de symétrie du maillage en fonction des modèles
que sur les efforts, bien qu’il ne soit pas trivial de comparer un calcul 2D à des mesures
3D. Les résultats obtenus sur la vitesse dans le rotor et le confinement servent à construire
le modèle de disque d’action empirique pour la Darrieus. Comme évoqué précédemment,
il n’existe pas de référence de comparaison pour le sillage de la turbine Darrieus dans le
sillage proche ou lointain. La simulation complète a pour but de servir de référence à la
construction du modèle simple. Dans la perspective d’introduire dans les modélisations
des valeurs de turbulence élevées couramment rencontrées en milieux naturels, on présente
ci-après une étude paramétrique des deux modèles (v 2 −f et kω −SST ) sur une large plage
de niveaux de turbulence amont. On étudiera successivement le cas complet rotor-stator
et le cas du modèle d’équivalence.
Le maillage du rotor dont le bon comportement vis-à-vis des détachements de couche limite
a été montré par les comparaisons expérimentales, est utilisé tout au long de cette partie
[33]. Plusieurs maillages pour la partie stator sont construits. L’objectif de cette étude
étant d’observer le rétablissement naturel de la vitesse dans le sillage, il est choisi de ne
pas agrandir les mailles près de la sortie pour ne pas dissiper artificiellement les structures.
2. Un second maillage étendu jusqu’à 15 diamètres est construit. Il est destiné à sur-
veiller l’évolution des structures cohérentes, dans le but de comprendre les phéno-
mènes d’instabilités observées. La taille de maille maximale d’environ ∆X = ∆Y =
D/40.
0.8 0.8
0.7 0.7
0.6 0.6
0.5 0.5
0.4 0.4
0.3 0.3
0.2 0.2
0.1 0.1
0 0
−0.1 −0.1
−0.2 −0.2
0 50 100 150 200 250 300 350 400 0 50 100 150 200 250 300 350 400
a) b)
1 1
0.9 0.9
0.8 0.8
0.7 0.7
0.6 0.6
0.5 0.5
0.4 0.4
0.3 0.3
0.2 0.2
0.1 0.1
0 0
−0.1 −0.1
−0.2 −0.2
−40 −30 −20 −10 0 10 20 30 40 −40 −30 −20 −10 0 10 20 30 40
c) d)
1 1
0.9 0.9
0.8 0.8
0.7 0.7
0.6 0.6
0.5 0.5
0.4 0.4
0.3 0.3
0.2 0.2
0.1 0.1
0 0
−0.1 −0.1
−0.2 −0.2
−40 −30 −20 −10 0 10 20 30 40 −40 −30 −20 −10 0 10 20 30 40
e) f)
Figure 4.11 – Colonne de gauche modèle kω − SST , droite v 2 − f . (a,b) CP sur une pale,
(c,d) CL projeté sur la tangente au cercle de rotation (force motrice), (e,f) CD projeté
(force résistante)
148 CHAPITRE 4. ETUDE DU SILLAGE DE LA TURBINE ACHARD
0.8 0.8
0.7 0.7
0.6 0.6
0.5 0.5
0.4 0.4
0.3 0.3
0.2 0.2
0.1 0.1
0 0
−0.1 −0.1
−0.2 −0.2
0 50 100 150 200 250 300 350 400 0 50 100 150 200 250 300 350 400
Figure 4.12 – Coefficient de puissance sur une pale, sur un tour, séparation des forces
totales (traits pleins) et des forces de pression (marqueurs). Gauche kω − SST , droite
v2 − f
la corde. Cette valeur faible entraîne suffisamment de dissipation à l’amont pour n’avoir
pratiquement pas de turbulence en entrée de turbine. Pour cette nouvelle étude le diamètre
hydraulique est modifié. On utilise la nouvelle valeur DH = 0, 36m qui correspond au
diamètre hydraulique de la section de la veine d’essai. L’intensité turbulente prend les
valeurs It = 0, 8; 8; 25%.
La puissance moyenne sur les trois pales en fonction du modèle et des paramètres d’en-
trée est reportée dans le tableau 4.1. On remarque tout d’abord que pour un diamètre
hydraulique faible, la valeur de k en entrée du calcul a peu d’influence sur les résultats. On
constate en effet une chute rapide de k entre l’entrée et la zone du rotor. Les pales voient
alors un écoulement avec une viscosité turbulente qui varie peu avec It . L’effet est plus
important pour une valeur de DH = 0, 36m. On constate un évolution opposée selon les
deux modèles : le CP moyen diminue fortement en v 2 −f alors qu’il augmente en kω −SST .
Afin de comprendre la différence observée sur la puissance, on s’intéresse à la puissance sur
une pale, mais également aux coefficients de portance et de trainée en fonction de l’incidence
de la pale. Les figures 4.11.a) et b) présentent le coefficient de puissance récupérée sur une
pale en fonction de la position angulaire. En modèle kω − SST , l’influence de l’intensité
turbulente en entrée est vue pour Dh = 0.36m à partir de It = 8%. Elle provoque une
augmentation du pic principal, en retardant la chute du couple. Le couple est le même
dans la zone θ = 0 − 80◦ . Dans la partie descendante, il est augmenté de 10 − 20% pour
It = 25% et la valeur minimum est plus élevée. L’augmentation est également vue dans la
deuxième zone motrice mais de façon moindre. En modèle v 2 −f , l’effet de la turbulence est
4.4. ETUDE NUMÉRIQUE DU SILLAGE EN MAILLAGE ROTATIF 149
CP 3 pales IT (%)
modèle DH (m) 0,8 8 25
v2 − f 0,032 0,53 0,47 0,42
0,36 0,42 0,27 0,16
kω − SST 0,032 0,31 0,32 0,35
0,36 0,31 0,36 0,47
Table 4.1 – Coefficient de puissance moyen pour la turbine, calculé sur un tour après 20
tours de convergence. On remarque l’influence opposée de l’intensité turbulente pour les
deux modèles kω − SST et v 2 − f .
très différent. L’amplitude du premier pic est drastiquement diminuée sur l’ensemble de la
rotation. Cette diminution se constate également à DH = 0, 032. Les figures 4.11 c,d) et e,f)
présentent respectivement les coefficients de portance et de trainée, projetés sur la tangente
au cercle de rotation. Ces grandeurs représentent la contribution motrice ou résistante pour
la machine. Concernant les coefficients de portance, l’influence de It est faible. Elle montre
une légère diminution de la portance maximale dans la zone motrice amont pour les deux
modèles. Une amélioration d’un facteur 3 du CL est constatée en modèle kω − SST dans
la deuxième phase du fait qu’on ne constate pas de deuxième décrochage.
L’effet de l’augmentation de l’intensité turbulente est surtout remarqué sur le coefficient
de trainée résistante. En modèle kω − SST , le calcul à faible intensité turbulente montre
une augmentation brusque du coefficient de trainée à α ≈ 25◦ pour It = 0, 8% et 8%,
qui n’apparaît pas à It = 25%. Le modèle v 2 − f montre une augmentation globale du
coefficient de trainée. La différence est attribuable en partie à la répartition du couple
entre les forces de pression, des contraintes visqueuses et turbulentes. On peut voir sur la
figure 4.12 le couple sur une pale dû aux forces de pression et la contribution totale. On
remarque que les efforts de trainée visqueuse sont négligeables en modèle kω − SST , alors
qu’elle impose une diminution significative du couple en modèle v 2 − f . De plus, la part de
la trainée visqueuse augmente avec l’intensité turbulente.
Le champ de vorticité autour de la pale est représenté pour les deux modèle de turbulence
pour deux positions angulaires sur la figure 4.13. La dynamique du décrochage en fonction
de It est très similaire en modèle v 2 − f , avec une intensité de la circulation plus forte à
It = 0, 8% et un décrochage qui semble intervenir légèrement plus tôt. En modèle kω −
SST , l’augmentation de la turbulence fait que le décrochage dynamique est retardé. Le
tourbillon de bord d’attaque reste accroché plus longtemps au profil, ce qui prolonge la
phase motrice. De plus, l’intensité du tourbillon est beaucoup plus forte avec le modèle
kω − SST . Néanmoins, dans ce cas l’écoulement montre la formation d’un tourbillon en
150 CHAPITRE 4. ETUDE DU SILLAGE DE LA TURBINE ACHARD
sens contraire entre le tourbillon principal et le bord d’attaque. Ce phénomène n’est pas
rapporté dans les différentes observations expérimentales.
θ = 110◦
It = 0, 8%
It = 25%
θ = 150◦
It = 0, 8%
It = 25%
kω − SST v2 − f
Figure 4.13 – Champ de vorticité autour de la pale, colonne de gauche, modèle kω −SST ,
droite, v 2 − f
152 CHAPITRE 4. ETUDE DU SILLAGE DE LA TURBINE ACHARD
Figure 4.14 – Contour des isovaleurs du taux de rotation dans les sillages des simulations
en modèle v 2 − f (haut) et kω − SST (bas)
place dans un domaine relativement confiné (L/D = 4). A titre d’exemple, un maillage de
largeur L/D = 20 est construit pour simuler une hydrolienne en milieu quasi-infini. On
choisit le modèle kω − SST qui permet de faire converger les calculs plus rapidement avec
une intensité d’entrée de 0, 8%. Le champ de vitesse transverse est tracé sur la figure 4.16.
On observe des fluctuations plus fortes que pour les mêmes paramètres en domaine confiné,
ainsi qu’un élargissement du sillage.
L’instabilité augmente en intensité quand l’énergie turbulente ambiante diminue. Une ex-
plication probable à cet effet est la nature instable des couches de mélange. Les faibles
valeurs de l’énergie cinétique turbulente limitent la diffusion turbulente, ce qui prolonge
l’existence de cet équilibre entre la zone déficitaire, et l’écoulement contournant. Pour le
cas de la figure 4.15.d, on observe des recirculations d’un ordre de grandeur du diamètre de
la turbine. Dans les autres cas, aucune recirculation n’est constatée mais uniquement un
mouvement transverse. Pour la simulation v 2 − f , le phénomène ressemble à une instabilité
4.4. ETUDE NUMÉRIQUE DU SILLAGE EN MAILLAGE ROTATIF 153
a)
b)
c)
d)
Figure 4.15 – Le champ de vitesse transverse instantanée dans le sillage est révélateur
de l’instabilité de l’écoulement qui diminue avec l’intensité turbulente ambiante. a,b,c)
kω − SST , respectivement It = 0, 8; 8; 25%, d) v 2 − f , It = 0, 8%
154 CHAPITRE 4. ETUDE DU SILLAGE DE LA TURBINE ACHARD
Figure 4.16 – Champ de vitesse transverse de la turbine en milieu infini pour It = 0.8%
avec le modèle de turbulence kω − SST
1.4 0.18
1.3
0.16
1.2
1.1 0.14
1
0.12
0.9
0.8 0.1
0.7
0.6 0.08
0.5
0.06
0.4
0.3 0.04
0.2
0.02
0.1
0 0
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1
a) b)
2 0.06
1.8 0.055
0.05
1.6
0.045
1.4
0.04
1.2 0.035
1 0.03
0.8 0.025
0.02
0.6
0.015
0.4
0.01
0.2 0.005
0 0
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1
c) d)
Figure 4.17 – Sillage proche comparés pour les modèles kω − SST et v 2 − f . La produc-
tion turbulente dans le sillage proche est largement différente entre les deux modèles de
turbulence v 2 −f et kω −SST . De plus, on remarque que le kω −SST dissipe les structures
cohérentes (c), illustrée par la variance du champ moyen (d), de façon plus rapide.
156 CHAPITRE 4. ETUDE DU SILLAGE DE LA TURBINE ACHARD
par la formule (4.18). Cette grandeur augmente la viscosité apparente de l’écoulement. Elle
représente un coefficient de diffusion de la quantité de mouvement.
k2
µt = C µ (4.18)
ǫ
La seconde est révélatrice des grandes structures déterministes qui provoquent une pertur-
bation du champ moyen. Dans notre cas, elle est importante dans les zones de transport
des vortex et lors du battement du tube de courant. La comparaison des deux modèles
montre un comportement opposé par rapport à ces grandeurs. Comme évoqué plus tôt le
modèle kω − SST dissipe les structures cohérentes rapidement, leur énergie se retrouve
dans l’énergie cinétique turbulente k (figure 4.17.a). Le niveau de k continue d’augmenter
pour être multiplié par 4 à X/D = 2, 5. Cette augmentation est attribuable à la fin de
la dissipation des tourbillons et aux forts gradients entre la partie déficitaire du sillage et
l’extérieur. La viscosité turbulente augmente de la même façon. Elle atteint le même ordre
de grandeurs que la viscosité dynamique. En modèle v 2 − f , la partie résolue de l’énergie
décroit jusqu’à une valeur comparable à celle de k après 2,5 diamètres. On peut voir sur
la figure 4.18 la valeur minimale du champ de vitesse moyen, en fonction de la distance au
rotor, pour différentes intensités turbulentes. On remarque que le calcul v 2 − f ne présente
pas de récupération de la vitesse dans un premier temps. Du fait de la faible viscosité
apparente, il y a peu de diffusion de la quantité de mouvement. Le fort gradient de vitesse
latéral se conserve donc sur quelques diamètres ce qui déclenche l’instabilité. Elle est res-
ponsable d’un mélange rapide dans la zone X/D = 6 − 10. Après 10 diamètres la vitesse
minimum montre un palier à ≈ 75% de la vitesse amont. L’augmentation du minimum de
la vitesse moyenne est attribuable au mélange par les grandes structures et non pas à une
diffusion de quantité de mouvement. La sur-dissipation du modèle v 2 − f dans le sillage
impose un sillage non-diffusif.
On rappelle que les calculs montrent une faible influence de l’intensité turbulente sur les
efforts pour DH = 0, 032 en raison de la dissipation rapide de k avant le rotor. Dans le
sillage, cette influence est également faible comme le montre la figure 4.18. Pour cette
raison, on réduira par la suite notre analyse au seul paramètre sur l’intensité turbulente
en prenant DH = 0, 36. Les figures 4.19 présentent la vitesse minimum (a) et la valeur
moyenne de k, K et µt dans le sillage de l’hydrolienne pour les trois intensités turbulentes
en modèle kω−SST . Elles servent à analyser et comprendre les différences dans la longueur
de réalimentation.
Pour It = 0, 8%, le niveau de turbulence est beaucoup plus bas à la sortie du rotor. La
différence d’intensité turbulente n’est pas compensée par la présence de la turbine. Le
niveau de turbulence atteinte à 15D est donc entièrement du à la production dans les zones
cisaillées.
4.4. ETUDE NUMÉRIQUE DU SILLAGE EN MAILLAGE ROTATIF 157
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30
1 4
0.9 3.5
0.8
3
0.7
2.5
0.6
2
0.5
1.5
0.4
1
0.3
0.2 0.5
0.1 0
0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30
a) b)
0.6 8
0.55
7
0.5
0.45 6
0.4
5
0.35
0.3 4
0.25
3
0.2
0.15 2
0.1
1
0.05
0 0
0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30
c) d)
Figure 4.19 – Vitesse relative et énergie cinétique turbulente modélisée dans le sillage de
la Darrieus simulée en modèle kω − SST - influence de la turbulence l’inlet
4.5. ETUDE DU SILLAGE PAR LES MÉTHODES DES TERMES SOURCES 159
à l’esprit que les simulations sont en 2D. En 3 dimensions, la longueur du sillage sera
plus courte, en raison de la diffusion de quantité de mouvement par les frontières
supérieures et inférieures.
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
−0.1
−0.2
0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30
X
impact a cette nouvelle zone de production. On impose une nouvelle forme à la distribution
de charge, qui permet une convergence plus rapide du modèle d’équivalence. La nouvelle
fonction de forme conserve les deux zones de travail amont et aval, mais impose maintenant
une force de trainée constante dans la direction Y, non nulle pour −0, 3 < ycell /D < 0, 3.
L’augmentation de la turbulence a un effet direct sur le fonctionnement du modèle. La
diffusion accélérée de quantité de mouvement provoque une réalimentation précoce de la
vitesse interne, comme on peut le voir sur la figure 4.22. L’effet de blocage est donc diminué,
ce qui améliore le CP donné par le modèle itératif. Le modèle donne une puissance et une
trainée qui augmente avec l’intensité turbulente. Les résultats obtenus en modèle kω −SST
sont donnés ci-dessous ; pour respectivement It = [0, 8; 8; 25] :
La vitesse minimum dans le sillage en fonction de la distance au rotor est tracée sur
la figure 4.20. On constate ici encore une augmentation de la vitesse de réalimentation
lorsque la turbulence augmente. Néanmoins cette influence est moins marquée que dans
la situation sans turbine pour lequel le cisaillement est plus faible. La différence entre les
modèles est également non négligeable, avec un écart de 10% sur la vitesse axiale à 15D.
On constate une différence beaucoup moins importante entre les différents niveaux de
turbulence que dans le cas cisaillé simple. On a vu dans le cas du calcul sans turbine que
4.5. ETUDE DU SILLAGE PAR LES MÉTHODES DES TERMES SOURCES 161
8e−02
7e−02
6e−02
5e−02
4e−02
3e−02
2e−02
1e−02
0e00
0 5 10 15 20 25 30
Figure 4.21 – Energie cinétique turbulente k moyenne sur la veine dans le sillage du
modèle de turbine
la dissipation de k est la même pour les différents modèles jusqu’à 4 diamètres. Le niveau
de turbulence au niveau de la turbine est donc le même entre les modèles de turbulence.
Les résultats intéressants que l’on peut tirer des profils présentés sur la figure 4.21 sont les
suivants : premièrement, le modèle kω − SST produit plus d’énergie cinétique turbulente
que les deux autres modèles, conçus pour limiter la production trop importante du modèle
k − ǫ de base. L’augmentation de l’énergie cinétique turbulente commence dans le rotor. Ce
qui n’est pas montré par la simulation complète ou les mesures LDV. Le modèle k − ǫ-PL
limite la production dans le zone très cisaillée des termes sources mais dissipe cette énergie
moins vite que le modèle v 2 −f . Une augmentation de la turbulence ambiante compense les
différences. La vitesse de réalimentation est en effet la même à It = 25%. Troisièmement,
le taux de turbulence après 25 diamètres est le même dans les 9 simulations.
Figure 4.22 – Champ de vitesse dans la zone des termes sources en modèle kω − SST en
fonction de l’intensité turbulente
une répartition en deux zones, permettant d’obtenir une vitesse moyenne dans la machine
correcte pour toute la plage de vitesse de rotation. L’influence de la répartition des forces
est mise en évidence dans cette partie. On compare trois répartitions :
2. Une seconde répartition en deux zones qui impose une force constante dans la direc-
tion transverse. Le terme source est nul pour |Y | > 5R/8 (rep.Const.Y)
3. Une répartition constante dans les deux directions. Le terme source est constant dans
la zone [|X| < R, |Y | < 5R/8] (rep.Uniforme).
Ces répartitions permettent d’obtenir une puissance et une force de trainée identiques avec
le schéma itératif de la méthode, pour un modèle de turbulence donné. Les trois précédentes
solutions sont utilisées pour simuler la turbine avec le modèle de turbulence v 2 −f . La vitesse
de réalimentation est très similaire avec les répartition 2 et 3, comme on peut le voir sur la
figure 4.23.a. La production turbulente plus importante dans ces deux cas stabilise le calcul
(figure 4.23.b), par contre le niveau de fluctuation résolue est très faible. La répartition 1
entraîne une instationnarité plus importante du sillage, qui apparaît dès l’aval immédiat
de la turbine (figure 4.24). Le démarrage de récupération du sillage suit les cas 1 et 2, mais
4.5. ETUDE DU SILLAGE PAR LES MÉTHODES DES TERMES SOURCES 163
0.9 0.04
0.8
0.035
0.7
0.03
0.6
0.025
0.5
0.4 0.02
0.3
0.015
0.2
0.01
0.1
0.005
0
−0.1 0
0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30
a) b)
0.016
0.014
0.012
0.01
0.008
0.006
0.004
0.002
0
0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30
1 0.4
0.9
0.35
0.8
0.7 0.3
0.6
0.25
0.5
0.4 0.2
0.3
0.15
0.2
0.1 0.1
0
0.05
−0.1
−0.2 0
0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30
Figure 4.25 – Effet de l’ajout des termes sources d’énergie turbulente et de dissipation
sur le sillage proposés par Rhétoré et al. [29]
4.6. COMPARAISON DES SILLAGES PROCHES 165
1.4 1.4
1.3 1.3
1.2 1.2
1.1 1.1
1 1
0.9 0.9
0.8 0.8
0.7 0.7
0.6 0.6
0.5 0.5
0.4 0.4
0.3 0.3
0.2 0.2
0.1 0.1
0 0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2
(a) (b)
0.6 0.6
0.55 0.55
0.5 0.5
0.45 0.45
0.4 0.4
0.35 0.35
0.3 0.3
0.25 0.25
0.2 0.2
0.15 0.15
0.1 0.1
0.05 0.05
0 0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2
(c) (d)
0.05 0.05
0.045 0.045
0.04 0.04
0.035 0.035
0.03 0.03
0.025 0.025
0.02 0.02
0.015 0.015
0.01 0.01
0.005 0.005
0 0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2
(e) (f)
Figure 4.26 – Profil de vitesse et d’énergie cinétique turbulente dans le sillage du modèle
de turbine
4.6. COMPARAISON DES SILLAGES PROCHES 167
nécessitent une validation par rapport à des données expérimentales. Une modification du
terme de production dans l’écoulement global est plus prédictive.
Points Clés
• La littérature concernant les sillages d’hydroliennes à axe vertical est peu fournie. En
particulier, il n’existe aucune référence expérimentale précise sur les turbine à flux
transverse.
• Les références à l’utilisation des modèles RANS avec le modèle du disque d’action
sont nombreuses. Les différents résultats et les corrections apportées au modèle k − ǫ
pour bien prendre en compte les sillages ne montrent pas de solution claire à adopter.
• Le sillage proche du modèle de turbine est analysé à l’aide de mesures LDV dans un
plan en 2 dimensions. La décomposition de l’énergie par la moyenne de phase donne
un outil cohérent de comparaison avec les calculs RANS.
• Le sillage proche montre une forte asymétrie due à l’intensité des tourbillons lâchés
par les pales. La symétrie semble être récupérée vers X/D=2,5.
• Deux type d’instabilités du sillage apparaissent dans les simulations : pour le modèle
v 2 −f , une instabilité des zones cisaillées intervient en formant des grandes structures
tourbillonnaires. Pour les autres simulations, une instabilité de jet est remarquée. Son
amplitude augmente avec la largeur du domaine disponible à l’écoulement.
• La réalimentation trop rapide pour les modèles de disques d’action vue dans le lit-
térature est due à l’intensité du cisaillement dans la zone des pertes de charge. La
réponse classique est d’imposer une augmentation de la dissipation artificielle, dans
une zone arbitraire ou avec des paramètres variables. Elle est peu satisfaisante car
peu prédictive.
Bibliographie
[1] A. Crespo, J. Hernandez, and S. Frandsen, “Survey of modelling methods for wind
turbine wakes and wind farms,” Wind energy, vol. 2, no. 1, pp. 1–24, 1999.
[3] L. Vermeer, J. N. Sørensen, and A. Crespo, “Wind turbine wake aerodynamics,” Pro-
gress in aerospace sciences, vol. 39, no. 6, pp. 467–510, 2003.
[4] B. Sanderse, S. Pijl, and B. Koren, “Review of computational fluid dynamics for wind
turbine wake aerodynamics,” Wind Energy, vol. 14, no. 7, pp. 799–819, 2011.
[5] G. I. Gretton, Hydrodynamic analysis of a vertical axis tidal current turbine. PhD
thesis, 2009.
[7] P. Alfredsson and J. Dahlberg, “A preliminary wind tunnel study of windmill wake
dispersion in various flow conditions,” Technical note AU-1499, part, vol. 7, 1979.
[8] V. P., “A wind tunnel study of the wake of a horizontal axis wind turbine,” 1978.
[9] A. Bahaj and L. Myers, “Shaping array design of marine current energy converters
through scaled experimental analysis,” Energy, vol. 59, pp. 83–94, 2013.
[12] D. Medici and P. Alfredsson, “Measurements on a wind turbine wake : 3d effects and
bluff body vortex shedding,” Wind Energy, vol. 9, no. 3, pp. 219–236, 2006.
170 CHAPITRE 4. ETUDE DU SILLAGE DE LA TURBINE ACHARD
[24] J. Mann, “The spatial structure of neutral atmospheric surface-layer turbulence,” Jour-
nal of Fluid Mechanics, vol. 273, pp. 141–168, 1994.
[29] P.-E. M. Rethore, N. N. Sørensen, A. Bechmann, and F. Zahle, “Study of the atmos-
pheric wake turbulence of a cfd actuator disc model,” 2009.
[30] A. El Kasmi and C. Masson, “An extended k–ε model for turbulent flow through
horizontal-axis wind turbines,” Journal of Wind Engineering and Industrial Aerody-
namics, vol. 96, no. 1, pp. 103–122, 2008.
[32] V. Aumelas, Modélisation des hydroliennes à axe vertical libres ou carénées : développe-
ment d’un moyen expérimental et d’un moyen numérique pour l’étude de la cavitation.
PhD thesis, Grenoble-INP, 2011.
[33] T. Maitre, E. Amet, and C. Pellone, “Modeling of the flow in a darrieus water turbine :
wall grid refinement analysis and comparison with experiments,” Jour of Renewable
Energy, accepted, 2012.
[34] E. Amet, “Simulation numérique d’une hydrolienne à axe vertical de type darrieus,”
2008.
172 CHAPITRE 4. ETUDE DU SILLAGE DE LA TURBINE ACHARD
Chapitre 5
Le placement relatif des turbines dans un parc hydrolien est un problème complexe. Les
sites choisis dans les applications réelles sont généralement contraints par la logistique
d’installation.
Plusieurs dispositions génériques sont envisagées ou développées par les chercheurs ou in-
dustriels. La figure 5.1 donne une illustration des trois principales solutions relevées par
Khan et al. [1]. L’ancrage au sol ou le lestage par système gravitaire est actuellement choisi
pour les implantations des prototypes industriels de plus grande taille, et en particulier
pour les turbines à flux axial. On peut citer les machines de technologie Open-Hydro ou
de la société Alstom. La topologie et la nature du lit est un critère de placement qui pré-
vaut pour l’instant sur la disposition pour avoir la meilleure puissance. Le montage sur
des structures flottantes est également utilisé. C’est une possibilité avantageuse en termes
de facilité de mise en place et de maintenance, et tout spécialement pour les turbines à
flux transverse verticales, car elle permet de maintenir la génératrice émergée en dehors
Figure 5.1 – Options de placement des turbines dans un canal, d’après Khan et al. [1]
173
174 CHAPITRE 5. APPLICATION DU MODÈLE À L’ÉTUDE DE PARCS
de l’eau. De plus, le placement dans une configuration particulière du parc est facilité et
la turbine fonctionne dans la partie haute de l’écoulement qui a la plus grande vitesse.
Néanmoins, comme l’illustre Khan, le partage d’usage impose de laisser une aire réservée
à la navigation. La fixation en bordure facilite largement l’installation et la maintenance
en solution, la production sera néanmoins sujette aux variations de niveau d’eau, et les
options de placement réduites aux zones de bord, ou les variations de vitesse et de hauteur
seront les plus importantes.
Tout en étant conscient de la complexité d’installation, il est primordial de disposer d’outils
pour l’optimisation des géométries de parc. Le modèle d’équivalence, dont la construction
a été décrite dans les précédentes parties de ce document a pour objectif de comprendre
l’interaction entre les différentes turbines, et son influence sur le rendement global du parc.
Différents effets vont entrer en jeu. De façon qualitative, on peut conjecturer les hypothèses
suivantes :
◦ Sur plusieurs lignes, le déficit de vitesse qui apparaît dans le sillage des turbines
amont a pour effet de diminuer l’énergie disponible pour les suivantes.
◦ L’intensité turbulente est augmentée par la forte production tourbillonnaire sur les
structures et le cisaillement entre les sillages et les zones de contournement.
Le dernier effet ne pourra être étudié avec le type de calcul effectué ici, dans lequel l’écou-
lement est forcé en vitesse. Pour les autres effets, plusieurs configurations théoriques de
parc sont étudiées avec l’objectif d’en tirer des enseignements généraux. L’ajout de cer-
taines compétences est envisageable dans le modèle afin de le rendre plus pertinent pour
l’optimisation de parc en situations réalistes. Certaines de ces améliorations sont proposées
dans un deuxième temps comme feuille de route de développements futurs.
5.1. RÉSULTATS GÉNÉRAUX 175
Symetry
Refined zone
Turbines
Inlet Outlet
180D
40D
6D
40D
Symetry
Nt
P
Pi
i
CP G = (5.1)
0.5ρU03 Nt A
ǫl = ∆A/D (5.2)
176 CHAPITRE 5. APPLICATION DU MODÈLE À L’ÉTUDE DE PARCS
0.54 0.75
0.52 0.7
0.5
0.65
0.48
0.6
0.46
0.44 0.55
0.42 0.5
0.4
0.45
0.38
0.4
0.36
0.34 0.35
0.32
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 2 4 6 8 10 12 14 16 18
a) b)
Figure 5.3 – CP global sur le parc pour une ligne de Nt machines à espacement et Nt
variable
0.75
0.7
0.65
0.6
0.55
0.5
0.45
0.4
0 1 2 3 4 5 6 7 8
Figure 5.4 – CP moyen par turbine en fonction de la distance à l’axe de symétrie du parc
Les dimensions de la figure de droite (Nt = 17) sont représentées avec un ratio de 5/17.
Cela permet d’évaluer la zone d’écoulement perturbé par rapport à la taille du parc. L’aug-
mentation du Reynolds de parc favorise la puissance. On constate de plus que la limite de
perturbation s’arrête assez loin des extrémités du domaine (non-visible ici), ce qui valide
l’hypothèse du milieu infini.
L’étude est menée sur des tailles de machines égales à celles du modèle réduit (R = 0, 0875
et U0 = 2, 3). Étant donné l’effet du nombre de Reynolds sur l’effet de blocage, il conviendra
d’effectuer des simulations sur des parcs aux dimensions non réduites
L’étude est mené sur des tailles de machines égale à celles du modèle réduit (R = 0, 0875 et
U0 = 2, 3). Étant donné l’effet du nombre de Reynolds sur l’effet de blocage, il conviendra
d’effectuer des simulations sur des parcs aux dimensions non réduites.
a) b)
Figure 5.5 – Carte de vitesse pour une ligne de 5 turbines (gauche) et 17 turbines. La
figure de droite est représentée avec une diminution de 5/17e
5.1. RÉSULTATS GÉNÉRAUX 179
tation de la trainée avec le confinement local sur une rangée de deux machines, et met en
avant que l’accélération créée entre les deux disques peut être récupérée avec un troisième
disque. D’après ces résultats, la présence de cet obstacle supplémentaire n’influence pas les
deux premiers. Cette configuration en quinconce (dite staggered en anglais) est également
trouvée comme la meilleure pour le rendement global par Gebreslassie [5] et Bai [6] dans
des études de parcs à grand nombre de machines.
Cette configuration est testée avec deux lignes consécutives : une première rangée de trois
turbines avec un espacement de ∆ Axe/D = 4, puis une seconde de deux turbines interca-
lées entre les premières. La distance longitudinale de la deuxième rangée (X2 ) est variable
afin de déterminer le point de meilleur rendement. On peut voir le résultat en terme de CP G
sur la figure 5.7. La situation dans laquelle les 5 machines sont alignées est trouvée comme
la meilleure. Le rendement baisse lorsque X2 augmente ou diminue. On peut comparer sur
la figure les valeurs obtenues avec le rendement global pour une ligne de 2 et 3 machines
isolées à espacement équivalent. Dans le cas infini (ie X2 très grand), on se rapproche de la
situation avec des parcs de respectivement 2 et 3 turbines en ligne. A partir de X2 /D = 4,
on passe sous ces limites. Les vues du champ de vitesse de la figure 5.6 montre qu’à partir
de cette distance la deuxième rangées reçoit directement l’impact du sillage de la première.
On remarque aussi que le placement latéral sur la deuxième ne semble pas optimal étant
donné l’élargissement du sillage, et qu’un décalage peut être intéressant si la distance entre
les lignes est imposée. Les configurations pyramidales convexe et concave de la figure sont
également testées. Le CP G ne montre pas d’augmentation par rapport au cas quinconce.
On explique ces résultats par le fait que la porosité du parc augmente pour une taille
caractéristique équivalente.
a) b)
c) d)
0.52
0.5
0.48
0.46
0.44
0.42
0.4
0.38
0.36
0.34
−2 −1 0 1 2 3 4 5 6 7 8 9 10 11 12
0.6
0.55
0.5
0.45
0.4
0.35
0.3
0.25
2 4 6 8 10 12 14 16 18 20 22
Figure 5.8 – CP moyen sur une ligne en fonction de la distance entre deux lignes
après la ligne 2. Cet effet est dû l’augmentation de l’énergie turbulente vue par la ligne 2.
Un effet de la turbulence sur le modèle, remarqué au chapitre précédent, est l’augmentation
du CP par une diffusion de quantité de mouvement rapide entre les deux zones de trainée.
Cet effet sur la puissance de l’intensité turbulente mérite d’être évalué et paramétré dans
le modèle pour pouvoir effectuer des calculs dans les sillages. L’inflexion de la courbe de
CP 1 vers 10 diamètres, est sans doute attribuable au fait que l’effet de l’énergie turbulente
diminue plus rapidement que le déficit de vitesse. On peut voir sur les figures 5.9.c) et
5.10.c), que l’intensité turbulente revient au niveau ambiant à cette distance, mais que
l’impact du sillage est toujours fort.
182
Figure 5.9 – Champ de vitesse dans le parc de trois lignes, la veine est de largeur Y/D=8, la distance
variable
CHAPITRE 5. APPLICATION DU MODÈLE À L’ÉTUDE DE PARCS
5.1. RÉSULTATS GÉNÉRAUX
➜ L’impact important de la turbulence sur le modèle, en créant une diffusion plus rapide
de la quantité de mouvement, et ainsi une augmentation du débit. La modélisation
U-RANS n’a pas permis d’avoir une idée claire sur l’effet de l’intensité turbulente sur
la performance. Il est donc indispensable de rentrer l’intensité turbulente comme pa-
ramètre du modèle, à partir de simulation plus précise (LES) ou de mesures. L’étude
bibliographique menée au chapitre précédent a cependant montré que la relation entre
intensité turbulence et puissance est difficile à estimer. On gardera en mémoire que
si ce paramètre est souvent mesuré dans les différentes expériences, la longueur de
mélange associée a une importance égale si ce n’est équivalente. Une illustration est
donnée par Blackmore et al. [7] par des mesures de trainée sur des disques poreux.
Cette force admet un minimum en fonction de la longueur de mélange pour des tailles
plus grandes et plus petites. Les effets combinés de l’intensité turbulente et de la lon-
gueur de mélange sont difficilement appréciables, et souvent peu abordés dans les
tests sur modèles réduits.
5.2. PISTES D’AMÉLIORATION DU MODÈLE EN VUE DE L’OPTIMISATION 185
Points Clés
• La puissance globale augmente avec le nombre de machines sur une ligne, mais éga-
lement avec la diminution de la distance latérale.
• En configuration quinconce, l’éloignement de la deuxième ligne diminue le rendement
global de façon monotone.
• L’augmentation de l’intensité turbulente produite par le sillage d’une turbine parti-
cipe à une réalimentation rapide du sillage des machines suivantes. Une forte intensité
turbulente tend à améliorer le rendement global du parc pour un nombre de rangée
supérieur à 2.
• L’influence du niveau de turbulence sur les turbines à flux transverse est fondamen-
talement importante et doit être prise en compte par le modèle simplifié.
• Un inventaire des améliorations pouvant être apportées à la méthode pour une
meilleur applicabilité est proposé.
186 CHAPITRE 5. APPLICATION DU MODÈLE À L’ÉTUDE DE PARCS
Bibliographie
[1] M. Khan, G. Bhuyan, M. Iqbal, and J. Quaicoe, “Hydrokinetic energy conversion sys-
tems and assessment of horizontal and vertical axis turbines for river and tidal appli-
cations : A technology status review,” Applied Energy, vol. 86, no. 10, pp. 1823–1835,
2009.
[2] I. G. Bryden and S. J. Couch, “Me1-marine energy extraction : tidal resource analysis,”
Renewable Energy, vol. 31, no. 2, pp. 133–139, 2006.
[3] T. Nishino and R. H. Willden, “Two-scale dynamics of flow past a partial cross-stream
array of tidal turbines,” Journal of Fluid Mechanics, vol. 730, pp. 220–244, 2013.
[4] L. Myers and A. Bahaj, “An experimental investigation simulating flow effects in first
generation marine current energy converter arrays,” Renewable Energy, vol. 37, no. 1,
pp. 28 – 36, 2012.
[6] G. Bai, J. Li, P. Fan, and G. Li, “Numerical investigations of the effects of different
arrays on power extractions of horizontal axis tidal current turbines,” Renewable Energy,
vol. 53, pp. 180–186, 2013.
[9] J. Steger, “The chimera method of flow simulation,” in Workshop on applied CFD, Univ
of Tennessee Space Institute, 1991.
Conclusion et perspectives
Le travail présenté tout au long de ce mémoire a eu pour but la construction d’un modèle
de turbines hydro-cinétiques à flux transverse de type Darrieus pour la simulation de parcs.
Les phénomènes physiques régissant la puissance récupérable par un tel système ont été
analysés en détail. Ceci permet la détermination des paramètres fondamentaux qui doivent
être pris en compte dans une représentation simplifiée.
Cette étude s’appuie sur des outils numériques et expérimentaux. D’un point de vue nu-
mérique, la prise en main et la validation d’une méthode de maillage rotatif dans le code
libre Code_Saturne permet la simulation de la turbine Achard en deux dimensions, dans
diverses configurations. Elle permet également le calcul de plusieurs turbines en géométrie
complète. Du côté expérimental, la mise en place d’un dispositif de mesure par vélocimétrie
laser Doppler (LDV) donne accès aux grandeurs cinétiques et turbulentes dans le sillage
proche d’un modèle réduit de la turbine.
Une représentation par des termes de quantité de mouvement, au sein d’un calcul d’écou-
lement, permet le calcul des performances de la turbine sur une large plage de paramètres
d’avance. Elle induit la réduction du temps de calcul de plusieurs ordres de grandeurs par
rapport à la géométrie complète. Une distribution des efforts dans la machine adaptée
aux turbines à flux transverse permet une paramétrisation du coefficient de puissance en
fonction de la vitesse interne. Le fonctionnement relatif de plusieurs modèles de turbulence
et leur réaction au niveau de turbulence ambiant est analysé et permet des conclusions
générales sur les modèles adaptés.
La capacité du modèle est illustrée par le calcul de plusieurs situations de parcs, mettant
en lumière une amélioration du rendement avec l’augmentation du nombre de machines.
187
188 Conclusion et perspectives
L’importance du taux de turbulence est également montrée comme fondamentale pour les
parcs de seconde génération (ie à plusieurs rangées de turbines). L’apport des différentes
parties est détaillé ci-après.
Le chapitre 1 expose une revue générale des concepts physiques des hydroliennes à flux
transverse. Grâce aux nombreux projet de recherche des dernières années, la documentation
est bien fournie, aussi ce chapitre est une synthèse de l’état de l’art. Il présente également
une revue des programmes d’essais expérimentaux qui ont été menés au LEGI depuis une
dizaine d’années.
durant la thèse de Ervin Amet (2008), avec le modèle v 2 − f cette fois. Elle montre une
très bonne comparaison avec les mesures PIV de Jonathan Bossard (2012) sur la dyna-
mique tourbillonnaire. On constate également une bonne cohérence sur le calcul du couple
instantané, bien que la comparaison simulation 2D avec les mesures 3D restent soumise à
interprétation. Le modèle v 2 − f prédit un décollement de la couche limite plus proche des
mesures que le modèle kω − SST . Le décrochage dynamique intervient plus tôt ce qui se
traduit par une baisse du coefficient de puissance moyen.
Des développements effectués dans Code_Saturne permettent de prendre en compte plu-
sieurs turbines et d’effectuer des calculs de parcs en géométrie complète.
Le troisième chapitre présente une revue des modèles simplifiés d’hydroliennes. Elle
montre que ces modèles se répartissent en deux principales catégories. La première utilise
une valeur empirique du coefficient de puissance. La puissance de la machine est déterminée
grâce au facteur d’induction déduit de la théorie idéale de Betz (méthode du disque d’ac-
tion). Dans le second cas, l’efficacité de l’hydrolienne est estimée à chaque calcul à partir
des vitesses locales dans le rotor et des coefficients aérodynamiques des pales (méthode
VBM).
Dans le présent travail, il est choisi d’aborder la construction du modèle en utilisant pleine-
ment les nombreuses données disponibles sur la machine et en simplifiant le fonctionnement
et l’implémentation dans un calcul d’écoulement. Le calcul en géométrie complète montre
l’augmentation de la puissance récupérée avec le confinement, mais également un décalage
du point de meilleur rendement. Un jeu de paramètres adimensionnels est construit en utili-
∗
sant la vitesse moyenne sur le disque. Les coefficients locaux de forces CX,Y et de puissance
CP∗ sont reliés au paramètre d’avance local λ∗ , pour obtenir une relation quasi-universelle.
L’étude montre que la simulation des forces de trainée par des termes puits de quantité
de mouvement permet le calcul correct de la vitesse interne pour les différents paramètres
d’avance et de confinement. Une fonction de forme analytique, basée sur la répartition
spatiale des efforts dans la turbine réelle permet l’utilisation du modèle sur un maillage
quelconque.
190 Conclusion et perspectives
Le fonctionnement du modèle est validé par comparaison avec les résultats de calcul en
maillage rotatif pour une turbine seule dans un canal. Un bon accord est constaté. Néan-
moins, le modèle surestime la puissance pour les fortes vitesses de rotation dans les cas très
confinés. Pour le cas de deux turbines en milieu infini une bonne correspondance est éga-
lement trouvée pour différentes vitesses de rotation pour un paramètre d’avance inférieur
à 2,5.
diffusion turbulente ce qui provoque une forte instabilité numérique. Les mauvais résultats
avec le modèle v 2 −f sont probablement attribuable à la nature 2D du calcul qui ne permet
pas la déstabilisation des tourbillons pas des fluctuations de vitesse verticales.
En ce qui concerne le modèle simplifié, la surproduction turbulente remarquée dans la
littérature est constaté avec le modèle k − ǫ. L’étude montre qu’elle est due aux zones très
fortement cisaillées dans la zone des termes sources et donc à la différence fondamentale
entre des représentations complètes et grossières. Le modèle k − ǫ-Production linéaire est
une modification du modèle k − ǫ classique. Elle impose une dépendance linéaire de la
production turbulente avec le taux de déformation de l’écoulement. Ce modèle permet
d’obtenir le comportement le plus semblable aux mesures dans le sillage proche.