Titre: Méthodologie d'optimisation des profils de vitesse à l'entrée des

Title: aspirateurs de turbines hydrauliques

Auteur:
Xinlin Li
Author:

Date: 2016
Type: Mémoire ou thèse / Dissertation or Thesis
Référence: Li, X. (2016). Méthodologie d'optimisation des profils de vitesse à l'entrée des
aspirateurs de turbines hydrauliques [Mémoire de maîtrise, École Polytechnique
Citation: de Montréal]. PolyPublie. https://publications.polymtl.ca/2413/

Document en libre accès dans PolyPublie

Open Access document in PolyPublie

URL de PolyPublie:
https://publications.polymtl.ca/2413/
PolyPublie URL:
Directeurs de
recherche: François Guibault, & Jean-Yves Trépanier
Advisors:
Programme:
Génie mécanique
Program:

Ce fichier a été téléchargé à partir de PolyPublie, le dépôt institutionnel de Polytechnique Montréal

This file has been downloaded from PolyPublie, the institutional repository of Polytechnique Montréal
https://publications.polymtl.ca
 UNIVERSITÉ DE MONTRÉAL

MÉTHODOLOGIE D’OPTIMISATION DES PROFILS DE VITESSE A L’ENTRÉE

DES ASPIRATEURS DE TURBINES HYDRAULIQUES

XINLIN LI

DÉPARTEMENT DE GÉNIE MÉCANIQUE

ÉCOLE POLYTECHNIQUE DE MONTRÉAL

MÉMOIRE PRÉSENTÉ EN VUE DE L’OBTENTION

DU DIPLÔME DE MAÎTRISE ÈS SCIENCES APPLIQUÉES

(GÉNIE MÉCANIQUE)

DÉCEMBRE 2016

© Xinlin Li, 2016.

 UNIVERSITÉ DE MONTRÉAL

ÉCOLE POLYTECHNIQUE DE MONTRÉAL

Ce mémoire intitulé:

MÉTHODOLOGIE D’OPTIMISATION DES PROFILS DE VITESSE A L’ENTRÉE

DES ASPIRATEURS DE TURBINES HYDRAULIQUES

présentépar : LI Xinlin

en vue de l’obtention du diplôme de : Maîtrise ès sciences appliquées

a étédûment acceptépar le jury d’examen constitué de :

M. GOSSELIN Frédérick, Ph. D., président

M. GUIBAULT François, Ph. D., membre et directeur de recherche

M. TRÉPANIER Jean-Yves, Ph. D., membre et codirecteur de recherche

M. GAUTHIER Maxime, M. Ing., membre

 iii

DÉDICACE

À mes chers parents, Zhiwei Li et Muchan Lin, qui m'aiment, croient en moi, m'inspirent et m'ont
soutenu àchaque étape du chemin.
 iv

REMERCIEMENTS

Tout d’abord, je tiens à exprimer ma sincère gratitude à mon directeur, le Professeur François
Guibault, pour son soutien constant et son mentorat tout au long du processus de recherche, pour
sa patience, sa motivation et son immense savoir. Je remercie également mon co-directeur le
Professeur Jean-Yves Trépanier pour ses suggestions et conseils utiles et perspicaces durant mes
études. Je remercie aussi très sincèrement les membres de mon jury de mémoire, le Professeur
Frederick Gosselin et Monsieur Maxime Gauthier.

En plus de mes directeurs, j’aimerais remercier le Docteur Christophe Devals pour ses précieux
conseils dans la mise en place et la validation du modèle de simulation CFD avec CFX et
OpenFOAM. Mes sincères remerciements vont également au Docteur Christophe Tribes et au
Professeur Sébastien Le Digabel pour leurs instructions afin de m’aider à comprendre l’algorithme
MADS et àutiliser NOMAD. J’ai également énormément apprécié l’aide de Madame Ying Zhang
afin de configurer l’outil de génération de maillages pour les deux cas de simulation de l’aspirateur.
Je remercie mes camarades de laboratoire pour les discussions stimulantes et pour tout le plaisir
que nous avons eu au cours des trois dernières années. Enfin, un merci spécial àMonsieur Thi Vu
d’Andritz Hydro, qui est un expert en optimisation de l’aspirateur, pour les discussions inspirantes
avec lui.
 v

RÉSUMÉ

La motivation àla base de ce projet découle des besoins exprimés par les ingénieurs-concepteurs
pour des outils spécialisés, répondant aux exigences des projets de réhabilitation de centrales
hydroélectriques. Dans un contexte de réhabilitation d’une centrale hydroélectrique existante, pour
concevoir une nouvelle roue de turbine afin de remplacer une ancienne roue et améliorer l'efficacité
énergétique globale de l'ensemble du système hydraulique, les concepteurs doivent déterminer quel
type d’écoulement en aval de la roue de turbine produira la moins grande perte d'énergie àl'intérieur
de l’aspirateur existant. L'approche proposée pour déterminer le comportement de l'écoulement
requis de l’aspirateur, consiste à formuler ce problème comme un problème d'optimisation des
conditions limites de vitesse à l’entrée de l’aspirateur.

Ce projet propose une méthodologie pour formuler et résoudre ce problème d'optimisation basée
sur l'algorithme d'optimisation MADS (Mesh Adaptive Direct Search) couplée àune approche de
simulation CFD basée sur la résolution des équations de Navier-Stokes en moyenne de Reynolds,et
utilisant le modèle de turbulence k-ε standard. Un cadre d'optimisation basésur Python, appelé
cfdOpt, a été développé pour mettre en œuvre cette méthodologie d'optimisation avec NOMAD et
OpenFOAM, qui sont des algorithmes d'optimisation et des codes de simulation CFD àsources
ouvertes. Une stratégie de parallélisation a été mise en œuvre dans cfdOpt pour profiter de la
capacitédes grappes de calcul àhautes performances afin d’accélérer le processus d'optimisation.

Pour garantir l'exactitude des simulations CFD, un exemple typique de simulation de l’aspirateur
basésur ANSYS CFX a étéutilisécomme référence pour configurer le cas de simulation avec
OpenFOAM. Le modèle de simulation CFD a également été validé par comparaison avec les
données expérimentales du projet Porjus U9.

La méthodologie a ététestée sur deux cas test distincts, un diffuseur conique et l’aspirateur Porjus
U9. Les résultats montrent que le facteur de perte d'énergie a étéréduit de plus de 60% dans les
deux cas d'optimisation par rapport au point de meilleur rendement (BEP) obtenu àl'aide du test
du corps solide en rotation. Ces résultats d'optimisation peuvent être utilisés comme référence de
conception pour les concepteurs de turbines travaillant sur des projets de réhabilitation de centrales
hydroélectriques.
 vi

ABSTRACT

The motivation at the root of this project stems from the need of design engineers for specialized
tools addressing the requirements of hydraulic power plant rehabilitation projects. To design a new
turbine runner to replace an old runner and improve the global energy efficiency of the whole
turbine system, designers must determine which types of downstream flow from the turbine runner
will yield the least energy loss inside the existing draft-tube. The proposed approach to determine
the required draft-tube flow behavior consists in formulating this as an inlet boundary condition
optimization problem.

This project proposes a methodology to formulate and solve this optimization problem based on
the Mesh Adaptive Direct Search (MADS) optimization algorithm coupled to an incompressible
Reynolds Averaged Navier-Stokes CFD simulation approach, using the standard k-ε turbulence
model. A Python-based optimization framework called cfdOpt was developed to implement this
optimization methodology with NOMAD and OpenFOAM, which are open-source optimization
algorithm and CFD simulation codes. A parallelization strategy was implemented in cfdOpt to take
advantage of high-performance cluster computing capacity for accelerating the optimization
process.

To guaranteed the correctness of the CFD simulations, a typical draft-tube simulation case based
on ANSYS CFX was used as a reference to setup the OpenFOAM simulation case. The CFD
simulation model was also validated by comparing with the experimental data from the Porjus U9
project.

The methodology was tested on two distinct test cases, a conical diffuser and the Porjus U9 draft-
tube. The results show that the energy loss factor was reduced by more than 60% in both
optimization cases compared with the best efficiency point found using the solid body rotation test.
These optimization results can be used as a design reference for turbine designers working on
rehabilitation projects of hydraulic power plants.
 vii

TABLE DES MATIÈRES

DÉDICACE ................................................................................................................................... III

REMERCIEMENTS ..................................................................................................................... IV

RÉSUMÉ........................................................................................................................................ V

ABSTRACT .................................................................................................................................. VI

TABLE DES MATIÈRES ...........................................................................................................VII

LISTE DES TABLEAUX .............................................................................................................. X

LISTE DES FIGURES .................................................................................................................. XI

LISTE DES SIGLES ET ABRÉVIATIONS .............................................................................. XV

LISTE DES ANNEXES ............................................................................................................. XVI

CHAPITRE 1 INTRODUCTION ............................................................................................... 1

CHAPITRE 2 REVUE DE LA LITTÉRATURE ....................................................................... 6

2.1. Paramétrage du profil de vitesse d’entrée de l’aspirateur ................................................ 9

2.1.1 Modèle analytique ...................................................................................................... 10

2.1.2 Modèle d’interpolation ............................................................................................... 11

2.2 Simulations CFD de l’écoulement interne dans un aspirateur ....................................... 14

2.2.1 Modèles de basse fidélité........................................................................................... 15

2.2.2 Modèles àhaute fidélité............................................................................................. 16

2.2.3 Traitement près de la paroi ......................................................................................... 22

2.3 Algorithmes d’optimisation numérique ......................................................................... 24

2.3.1 Méthodes basées sur les gradients .............................................................................. 24

2.3.2 Méthodes sans gradient .............................................................................................. 28

2.4 Résumé........................................................................................................................... 33

2.4.1 Paramétrage du profil de vitesse ................................................................................ 33

 viii

2.4.2 Simulation CFD.......................................................................................................... 34

2.4.3 Algorithme d’optimisation ......................................................................................... 35

CHAPITRE 3 MISE EN OEUVRE .......................................................................................... 38

3.1 Paramétrage .................................................................................................................... 39

3.1.1 Modèle analytique ...................................................................................................... 39

3.1.2 Modèle d’interpolation ............................................................................................... 39

3.2 Simulation CFD.............................................................................................................. 42

3.2.1 CFX ............................................................................................................................ 42

3.2.2 OpenFOAM ................................................................................................................ 42

3.2.3 Conditions aux limites ................................................................................................ 43

3.2.4 Traitement près de la paroi ......................................................................................... 45

3.2.5 Schéma de discrétisation ............................................................................................ 46

3.2.6 Initialisation du champ d’écoulement ........................................................................ 46

3.2.7 Évaluation de la fonction objectif .............................................................................. 47

3.2.8 Critères de convergence ............................................................................................. 50

3.2.9 Validation de la cohérence entre OpenFOAM et CFX .............................................. 50

3.3 Optimisation ................................................................................................................... 53

3.3.1 Point de départ de l’optimisation ............................................................................... 53

3.3.2 Contraintes non-linéaires............................................................................................ 54

3.3.3 NOMAD et parallélisation ......................................................................................... 55

CHAPITRE 4 RÉSULTATS .................................................................................................... 59

4.1 Cas du diffuseur conique ................................................................................................ 59

4.1.1 Maillage ...................................................................................................................... 60

4.1.2 Solution initiale pour l’optimisation .......................................................................... 66

 ix

4.1.3 Résultats de l’optimisation ......................................................................................... 67

4.2 Cas de l’aspirateur Porjus U9 ......................................................................................... 72

4.2.1 Maillage ...................................................................................................................... 72

4.2.2 Validation avec les données expérimentales .............................................................. 73

4.2.3 Solution initiale pour l’optimisation .......................................................................... 76

4.3 Résultats d’optimisation ................................................................................................. 78

4.3.1 Cas # 1 ........................................................................................................................ 78

4.3.2 Cas # 2 ........................................................................................................................ 82

4.3.3 Discussion .................................................................................................................. 86

CHAPITRE 5 CONCLUSION ET RECOMMANDATIONS ................................................. 87

5.1 Conclusion ...................................................................................................................... 87

5.2 Travail futur.................................................................................................................... 89

BIBLIOGRAPHIE ........................................................................................................................ 91

ANNEXES .................................................................................................................................... 96
 x

LISTE DES TABLEAUX

Tableau 3.1 Conditions aux limites ................................................................................................ 43

Tableau 4.1 Information concernant les maillages du diffuseur conique ...................................... 60

Tableau 4.2 Convergence de tous les cas de simulation dans la Figure 4.5 ................................... 65
 xi

LISTE DES FIGURES

Figure 1.1 Nombre cumulatif de barrages construits par région ...................................................... 1

Figure 1.2 Centrale hydroélectrique ................................................................................................. 2

Figure 2.1 Vortex forcé.................................................................................................................... 6

Figure 2.2 Organigramme d’optimisation de profile d’entrée de vitesse (Galvan, 2007) ............... 8

Figure 2.3 Système de coordonnées cylindriques pour la représentation des profils d’entrée ...... 10

Figure 2.4 Comparaison entre la spline cubique monotone et classique d’Hermite ...................... 13

Figure 2.5 Organigramme d’un algorithme génétique ................................................................... 30

Figure 2.6 L'étape sondage de la recherche de coordonnées.......................................................... 30

Figure 2.7 Recherche de coordonnées ............................................................................................ 31

Figure 2.8 Recherche de patrons généralisée ................................................................................. 32

Figure 2.9 Comparaison entre les directions de sondage de GPS (gauche) et MADS (droite) ..... 32

Figure 3.1 Diagramme de cfdOpt ................................................................................................... 38

Figure 3.2 Profil de vitesse typique pour les cas dtpmaker sans moyeu ........................................ 40

Figure 3.3 Profil de vitesse axiale corrigépour le débit................................................................. 41

Figure 3.4 Profil de vitesse tangentielle (dtpmaker) ...................................................................... 41

Figure 3.5 La partie d’extension du maillage de l’aspirateur ......................................................... 48

Figure 3.6 Le facteur de perte d’énergie (gauche) et les courbes résiduelles d’un cas typique de
CFD (droite) ........................................................................................................................... 50

Figure 3.7 Comparaison entre swak4foam et CFD-Post ................................................................ 51

Figure 3.8 Le profil de vitesse d’entrée pour la comparaison des conditions aux limites d’entrée52

Figure 3.9 Comparaison des profils d’entrée 𝑘 .............................................................................. 52

Figure 3.10 Comparaison des profils d’entrée 𝜀 ............................................................................ 53

Figure 3.11 Courbe résiduelle typique d’OpenFOAM................................................................... 55

 xii

Figure 3.12 Mode série NOMAD .................................................................................................. 56

Figure 3.13 Mode MPI NOMAD ................................................................................................... 57

Figure 3.14 Mode d’évaluation de bloc NOMAD ......................................................................... 58

Figure 4.1 Dimensions du diffuseur conique ................................................................................. 59

Figure 4.2 Maillage de diffuseur conique (Moyen) ....................................................................... 60

Figure 4.3 Courbes de facteur de perte d’énergie en fonction du nombre de tourbillon

(OpenFOAM) ......................................................................................................................... 61

Figure 4.4 Courbes de facteur de perte d’énergie en fonction du nombre de tourbillon (CFX) .... 61

Figure 4.5 Courbes de facteur de perte d’énergie et nombre de tourbillon (OpenFOAM et CFX)
................................................................................................................................................ 62

Figure 4.6 Comparaison entre les différents outils de post-traitement .......................................... 63

Figure 4.7 Courbes résiduelles typiques pour les cas de non-convergence dans CFX (gauche) et
OpenFOAM (droite) ............................................................................................................... 64

Figure 4.8 Courbes résiduelles typiques pour les cas de convergence CFX (gauche) et OpenFOAM
(droite) .................................................................................................................................... 64

Figure 4.9 Courbe du facteur de perte d’énergie en fonction du nombre de tourbillon (Diffuseur
conique) .................................................................................................................................. 66

Figure 4.10 Point de départ de l’optimisation ................................................................................ 67

Figure 4.11 Le graphique de convergence NOMAD (Diffuseur conique) .................................... 68

Figure 4.12 Solution optimale du diffuseur conique ...................................................................... 68

Figure 4.13 Courbe de facteur de perte d’énergie (Diffuseur conique) ......................................... 69

Figure 4.14 Courbes résiduelles (Diffuseur conique) .................................................................... 69

Figure 4.15 Champ de vitesse de la solution optimale (OpenFOAM) ........................................... 70

Figure 4.16 Contour de la pression de la solution optimale (OpenFOAM) ................................... 70

Figure 4.17 Champ de vitesse de la solution optimale (CFX) ....................................................... 71

 xiii

Figure 4.18 Contour de la pression de la solution optimale (CFX) ............................................... 71

Figure 4.19 Banc d’essai de l’aspirateur Porjus U9 (Mulu & Cervantes, 2009)............................ 72

Figure 4.20 Maillage de l’aspirateur Porjus U9 ............................................................................. 73

Figure 4.21 Profils de vitesse mesurés et numériques dans la section I......................................... 74

Figure 4.22 Profils de vitesse mesurés et numériques dans la section II ....................................... 74

Figure 4.23 Profils de vitesse mesurés et numériques dans la section III ...................................... 75

Figure 4.24 Contour de pression (OpenFOAM) ............................................................................ 75

Figure 4.25 Contour de pression (CFX) ......................................................................................... 75

Figure 4.26 Courbe de facteur de perte d’énergie en fonction du nombre de tourbillon (Porjus U9)
................................................................................................................................................ 76

Figure 4.27 Point de départ d’optimisation (Porjus U9) ................................................................ 77

Figure 4.28 Champ de vitesse du point de départ (Porjus U9) ...................................................... 77

Figure 4.29 Contour de pression du point de départ (Porjus U9) .................................................. 77

Figure 4.30 Graphique de convergence de l’optimisation (Cas # 1) .............................................. 78

Figure 4.31 Solution optimale (Cas # 1) ........................................................................................ 79

Figure 4.32 Courbes résiduelles de la solution optimale (Cas # 1) ................................................ 79

Figure 4.33 Courbe de facteur de perte d’énergie de la solution optimale (Cas # 1) ..................... 80

Figure 4.34 Champ de vitesse de la solution optimale (Cas # 1, OpenFOAM) ............................. 80

Figure 4.35 Champ de vitesse de la solution optimale (Cas # 1, CFX) ......................................... 81

Figure 4.36 Contour de pression la solution optimale (Cas # 1, OpenFOAM).............................. 81

Figure 4.37 Contour de pression la solution optimale (Cas # 1, CFX) .......................................... 81

Figure 4.38 Les gammes des points de contrôle tangentiels (Cas # 2) .......................................... 82

Figure 4.39 Graphique de convergence de l’optimisation (Cas # 2) .............................................. 83

Figure 4.40 Solution optimale (Cas # 2) ........................................................................................ 83

Figure 4.41 Courbes résiduelles de la solution optimale (Cas # 2) ................................................ 84

 xiv

Figure 4.42 Courbe de facteur de perte d’énergie de la solution optimale (Cas # 2) ..................... 84

Figure 4.43 Champ de vitesse de la solution optimale (Cas # 2, OpenFOAM) ............................. 85

Figure 4.44 Champ de vitesse de la solution optimale (Cas # 2, CFX) ......................................... 85

Figure 4.45 Contour de pression de la solution optimale (Cas # 2, OpenFOAM) ......................... 85

Figure 4.46 Contour de pression de la solution optimale (Cas # 2, CFX) ..................................... 86

 xv

LISTE DES SIGLES ET ABRÉVIATIONS

BEP Best Efficiency Point

CFD Computational Fluid Dynamics

CS Coordinate Search

GA Genetic Algorithm

GPS Generalized Pattern Search

MADS Mesh Adaptive Direct Search

RANS Reynolds-averaged Navier–Stokes

SBR Solid Body Rotation

 xvi

LISTE DES ANNEXES

Annexe A – L’EXPRESSION DU POST-TRAITEMENT ........................................................... 96

 1

CHAPITRE 1 INTRODUCTION

L’hydroélectricité est la plus ancienne source d’énergie renouvelable dans le monde. Les humains
ont exploité l’eau pour effectuer des travaux depuis des milliers d’années. Dans les premiers temps,
les gens ont construit leurs ateliers de travail àproximitéd’une rivière et ils y ont extrait de l'énergie
àtravers des roues hydrauliques. La première référence connue se trouve dans un poème grec de
l’an 85 avant Jésus-Christ. Cette ressource énergétique ancienne s’est développée rapidement en
Europe, en Amérique du Nord, et en Asie après la fin de la Seconde Guerre mondiale. La capacité
de production hydroélectrique installée dans le monde a augmenté de façon très significative,
passant de moins de 50 GW en 1950 àplus de 978 GW en 2009.

Figure 1.1 Nombre cumulatif de barrages construits par région1

En tant que générateur d’électricité, l’hydroélectricité est importante. À un niveau de 15,9% de la

production totale d’électricitédans le monde, l’hydroélectricité occoupe actuellement une place
plus importante que la puissance nucléaire, qui se situe à 13,5%. L’hydroélectricité est donc une
source d’énergie renouvelable très importante pour la production d’électricité.

En 2010, le Canada, avec une capacitéhydroélectrique installée de 74 GW, était le quatrième pays
au monde en termes de capacitéinstallée ("Global Hydro Power Report,") et près de la moitiéde

1
Source de données: (World Commission on, 2000)
 2

sa capacité est située dans la Province de Québec (36,5 GW). Ces centrales hydroélectriques
fournissent 98% de l’électricité au Québec, ce qui représente 40% de l’énergie totale consommée
dans cette province. Par conséquent, l’hydroélectricité joue un rôle important dans la société
québécoise.

Bien que la demande mondiale pour de l’énergie renouvelable soit en croissance continue et que
l’hydroélectricité joue un rôle vital dans cette croissance, le nombre de barrages construits
annuellement tend àdiminuer pour de multiples raisons, telles que la hausse du coût du travail et
la préoccupation croissante du public concernant l’impact environnemental des grands barrages.

Cette contradiction a incitéles ingénieurs àtrouver des moyens permettant d’extraire plus d’énergie
à partir de centrales hydroélectriques existantes au lieu d’en construire de nouvelles. La
réhabilitation des centrales hydroélectriques existantes - donc qui augmente la production d’énergie
tout en étendant leur durée de vie - devient de plus en plus importante. Le présent projet trouve sa
motivation dans les nouveaux besoins des ingénieurs pour des outils de conception adaptés àla
réhabilitation des centrales hydroélectriques.

Figure 1.2 Centrale hydroélectrique2

2
Source: File:Hydroelectric dam.svg - https://en.wikipedia.org
 3

Sans aucun doute, la turbine hydraulique est l’élément technologique central impliqué dans la
conversion de l’énergie hydraulique. Parmi les différents types de turbines hydrauliques, la turbine
Kaplan est la turbine axiale la plus largement utilisée dans le monde. Ce type de turbine permet
une grande production d’hydroélectricité de façon efficace dans le cas de conditions de débit élevé,
et de faibles hauteurs de chute.

L’aspirateur est l’un des composants les plus importants à l’intérieur de la turbine hydraulique. Il
convertit la pression dynamique de l'écoulement en pression statique par la décélération de
l'écoulement avant qu’il ne se déverse dans la rivière en aval. Il représente 20% à 50% de l’énergie
totale qui peut être récupérée dans une centrale électrique àfaible hauteur de chute (Andersson &
Dahlbäck, 1998). La performance de l’aspirateur dépend des distributions de vitesse à l’entrée de
la turbine et d’autres facteurs tels que la cavitation, la profondeur de l’eau en aval, le point de
fonctionnement de la turbine, la traînée, les décollements et les écoulements secondaires. Tous ces
facteurs dépendent non seulement de la forme géométrique de l’aspirateur mais ils dépendent aussi
fortement de la conception de la roue de la turbine.

Dans un projet de réhabilitation d’une centrale hydroélectrique, la bâche spirale et l’aspirateur sont
habituellement conservés en raison du fait qu’ils font partie du barrage et qu’ils sont généralement
construits en béton. Certains composants comme le générateur, les aubes directrices et la roue sont
remplacés. Par conséquent, échanger une roue nouvellement conçue qui correspond mieux à
l’aspirateur existant est le moyen le plus pratique et le plus efficace pour améliorer l’efficacité
globale de la production d’électricité de l’ensemble de la turbine.

Traditionnellement, pour obtenir la meilleure adéquation entre la roue et l’aspirateur, des essais sur
plusieurs modèles d’aspirateur et de roues étaient menées pour vérifier les designs des roues.
Toutefois, ces essais modèles sont très coûteux, de sorte que les concepteurs de turbines ne peuvent
pas explorer l’espace d’optimisation de façon approfondie et systématique. Ils peuvent seulement
ajuster certains paramètres de conception en fonction de leur propre expérience. Par conséquent, la
conception finale de la roue de turbine est généralement une conception pratique plutôt qu’une
conception optimale.

De nos jours, avec le développement rapide de l’informatique à haute performance et avec des
modèles CFD à haute-fidélité, les concepteurs sont en mesure d’obtenir à faible coût des
prédictions précise concernant l'écoulement dans l’aspirateur et de prévoir la performance de
 4

l’aspirateur sans faire des essais coûteuses. Une conception améliorée de la roue de turbine,
correspondant aux mieux àl’aspirateur existant, peut ainsi être obtenue sur la base de ces nouvelles
techniques. Des résultats précédents ont montré que le remplacement de la roue de la turbine,
obtenue grâce à des technologies modernes, tout en gardant les autres structures de turbines
existantes, permettait d’augmenter la production d’énergie d’une centrale hydroélectrique par un
facteur pouvant atteindre de 10 à30% (Galván, Reggio, & Guibault, 2012). Cependant, les études
actuelles sur l’optimisation des turbines hydrauliques sont encore principalement axées sur
l’optimisation de la forme de la pale de la turbine et l’optimisation géométrique de l’aspirateur,
dont seules les futures centrales hydroélectriques bénéficieront. Plus d’efforts devraient être voués
à l’amélioration de la conception de turbines hydrauliques pour des installations existantes.

Dans un projet de réhabilitation, la première étape consiste à déterminer quel type d’écoulement
en aval de la roue de turbine peut réduire la perte d’énergie dans l’aspirateur existant et maximiser
son efficacitéde récupération. Cette analyse est également connue comme l’optimisation du profil
de la vitesse d’entrée pour l’aspirateur. Les résultats de la résolution de ce problème d’optimisation
seront utilisés comme un objectif de conception de la nouvelle roue de turbine.

Dans une étude récente, Galvan a présenté une méthodologie d’optimisation pour les profils de
vitesse d’entrée des aspirateurs basée sur une spécification analytique des conditions limites de
vitesse à l’entrée d’un cône d’aspirateur (Galvan, 2007). Dans cet article, Galvan a développéune
méthodologie d’optimisation impliquant trois logiciels commerciaux, nommément iSIGHT,
MATLAB et FLUENT.

Le présent projet vise à améliorer la méthodologie d’optimisation proposée par Galvan, en y

intégrant d’une part une méthode plus souple de représentation des profils de vitesse afin d’élargir
l’espace d’optimisation et d’autre part en utilisant un algorithme d’optimisation numérique plus
efficace pour accélérer le processus. Au lieu d’un logiciel commercial, cette nouvelle méthodologie
d’optimisation sera mise en œuvre avec ses homologues sources ouvertes pour éviter des frais de
licences coûteuses. Par conséquent, le problème d’optimisation pourra être résolu àune plus grande
échelle et devrait être résolu plus rapidement que précédemment puisque le nombre maximal de
cas de simulation n’est plus limité par les licences de logiciels commerciaux.

Ce mémoire est organisécomme suit : le chapitre 2 présente la revue de la littérature concernant

l’optimisation de la vitesse d’entrée de l’aspirateur et justifie les décisions prises au cours de ce
 5

projet. Une méthode d’optimisation améliorée est identifiée. Le chapitre 3 est consacréàla mise
en œuvre de cette méthodologie. Un code python appelé cfdOpt a été développé pour intégrer deux
logiciels différents et construire la boucle d’optimisation. Les résultats de deux cas test pour des
aspirateurs différents sont présentés au chapitre 4, suivi par la conclusion de ce projet au chapitre 5.
 6

CHAPITRE 2 REVUE DE LA LITTÉRATURE

Dans une méthodologie d’optimisation traditionnelle d’une turbine hydraulique, chaque composant
est conçu de façon indépendante. Afin de tenir compte de la présence de l'aspirateur, lors de la
conception de la roue, une approche pour découpler les composants consiste àposer l’hypothèse
que le débit en aval de la roue de turbine - qui est aussi le débit d’entrée de l’aspirateur – est
approximativement décrit par une combinaison d’écoulement axial uniforme et de vortex forcé
(Figure 2.1), aussi connu comme un écoulement de type "corps solide en rotation".

Figure 2.1 Vortex forcé

L’écoulement tourbillonnaire du corps solide peut être décrit directement par sa vitesse axiale et sa
vitesse angulaire. Cependant, dans le domaine de la conception de turbines hydrauliques, les
ingénieurs préfèrent définir un écoulement d’entrée de l’aspirateur par son débit et le nombre de
tourbillon S (Swirl number), qui est un nombre sans dimension indiquant le rapport entre la quantité
de mouvement tangentielle et la quantitéde mouvement axiale de l’écoulement.

𝑅
∫0 (𝜌𝑉𝑎 )(𝑟𝑉𝑡 )𝑟 𝑑𝑟
𝑆= 𝑅 2-1
𝑅 ∫0 (𝜌𝑉𝑎 )(𝑟𝑉𝑎 ) 𝑑𝑟

Pour un point d'opération donnéde la turbine le débit d’entrée est généralement fixé. La pratique a
prouvéque la performance de l’aspirateur dépend en grande partie du nombre de tourbillon de
l’écoulement d’entrée. Pour une géométrie spécifique de l’aspirateur, il y a un point de rendement
optimal (Best Efficiency Point ou BEP) qui est fonction du nombre de tourbillon que la
méthodologie traditionnelle d’optimisation de l’aspirateur essaie de trouver (Figure 4.9).
 7

Dans le processus d’optimisation de la turbine, les concepteurs tentent d’extraire la majeure partie
de l’énergie cinétique de rotation de l’écoulement avec la roue de la turbine, de sorte qu’un
tourbillon de faible intensitéen aval de la roue est recherché. En pratique, ni un écoulement sans
tourbillon, ni un écoulement avec trop de tourbillon ne peuvent atteindre un rendement élevéà
l’intérieur de l’aspirateur. Les deux types d’écoulement réduisent la performance globale de la
turbine hydraulique, et le niveau approprié de tourbillon dans l’écoulement doit donc être
soigneusement identifié.

Les ordinateurs modernes et des techniques de simulation numérique permettent aux chercheurs de
mieux comprendre les problèmes complexes d’optimisation en ingénierie impliquant des
écoulements fluides. Grâce à l’élaboration de techniques CFD, le coût de la prédiction de la
performance de l’aspirateur a été considérablement réduit, le goulot de la méthode d’optimisation
n’est plus le nombre limité d’essais. Les concepteurs peuvent maintenant résoudre des problèmes
d’optimisation comportant un nombre important de variables. Par conséquent, plusieurs nouvelles
méthodologies d’optimisation ont étéproposé. Au cours des dernières années plusieurs études ont
présentédes méthodologies d’estimation et d’optimisation de la performance dans les aspirateurs.

(Marjavaara & Lundström, 2006) ont proposé une méthodologie d’optimisation qui permet
d’améliorer la géométrie du talon de l’aspirateur pour maximiser le facteur de récupération de
pression moyenne et minimiser le facteur de perte d’énergie.

(Puente, Reggio, & Guibault, 2003) ont décrit une stratégie d’optimisation automatique pour les
sections de l’aspirateur d’une turbine hydraulique en combinant CFX et iSIGHT.

Auparavant, la plupart des optimisations d’aspirateurs portaient sur l’optimisation de la forme de

l'enveloppe, qui ne peut être mise en œuvre que dans les futures centrales hydroélectriques puisque
les formes d’aspirateurs sont généralement non modifiables une fois que la centrale hydroélectrique
a été construite. Cependant, à partir de la Figure 1.1, on peut observer que la période de
développement rapide de barrages hydroélectriques semble révolue. Les projets de réhabilitation
des centrales hydroélectriques existantes, qui prolongent leur durée de vie et améliorent leur
performance, deviennent de plus en plus importants pour l’hydraulique industrielle.

Dans un projet de réhabilitation d’une turbine hydraulique, la bâche spirale et l’aspirateur sont
généralement conservés et d’autres composants comme le générateur, des aubes directrices et la
roue seront mis àjour avec des modèles améliorés. Pour concevoir ces nouveaux composants qui
 8

correspondent mieux à l’aspirateur existant, les informations sur l’écoulement à la section de

connexion entre la roue et l’aspirateur sont vitales. Cette information peut être obtenue en résolvant
le problème d’optimisation du profil de vitesse d’entrée pour l’aspirateur existant.

Dans l’optimisation de l’écoulement àl’entrée de l’aspirateur, de nos jours, les concepteurs de

l’aspirateur vont plus loin et mettent dorénavant l’emphase sur l’optimisation des profils de vitesse
complets plutôt que seulement sur le nombre de tourbillon. (Galvan, 2007) a proposé une
méthodologie d’optimisation pour résoudre le problème de profil d’optimisation de la vitesse
d’entrée pour l’aspirateur, qui comprend trois parties principales différentes, à savoir le
paramétrage du profil de vitesse, la simulation de l’écoulement et l’algorithme d’optimisation
(Figure 2.2).

Figure 2.2 Organigramme d’optimisation de profile d’entrée de vitesse (Galvan, 2007)

Au début de chaque boucle itérative, l’algorithme d’optimisation fournit un ensemble de variables

de conception. Ensuite, ces variables de conception sont envoyées àla partie de paramétrage du
profil de vitesse, oùles variables de conception sont converties en conditions aux limites de vitesse
àl’entrée. Les conditions aux limites d’entrée sont lues au début de la simulation et le code CFD
résout les équations de Navier-Stokes pour prédire le comportement de l’écoulement. En fin de
compte, les performances de cet ensemble de variables de conception sont retournées à l’algorithme
d’optimisation pour fermer la boucle d’optimisation.
 9

L’algorithme d’optimisation de la méthodologie de Galvan est la mise en œuvre de l’algorithme

génétique multi-îles fourni dans le logiciel commercial iSIGHT, développépar Dassault Systèmes.
Ce logiciel permet d’intégrer et de connecter différents codes pour la construction de boucles
itératives afin de résoudre des problèmes d’optimisation en ingénierie. Dans sa méthodologie,
Galvan utilise une méthode de paramétrage du profil de vitesse basée sur un modèle analytique qui
fut proposépar Susan-Resiga et al. (Romeo Susan-Resiga, Ciocan, Anton, & Avellan, 2006), Cette
approche sera décrite plus en détails àla section 2.1.1. Le facteur de perte d’énergie de l’aspirateur
a étésélectionnécomme fonction objectif pour l’optimisation. Ce facteur est défini comme,

1 1
𝐴𝑖𝑛 ∫𝑖𝑛 𝑃𝑡 𝑑𝐴 − 𝐴𝑜𝑢𝑡 ∫𝑜𝑢𝑡 𝑃𝑡 𝑑𝐴
𝜁= 2-2
1 𝑄 2
2 𝐴𝑖𝑛 )
𝜌 (

Où𝑃𝑡 est la pression totale, 𝐴𝑖𝑛 et 𝐴𝑜𝑢𝑡 sont les sections d’entrée et de sortie et 𝑄 est le débit
volumétrique.

En fait, le facteur de perte d’énergie est la différence entre la pression totale moyennée àla section
d’entrée et de sortie de l’aspirateur, normalisée par l’énergie cinématique àl’entrée de l’aspirateur.
Il mesure le pourcentage total de perte de pression à l’intérieur de l’aspirateur que les concepteurs
d’aspirateur tentent de réduire.

Outre Galvan, (R Susan-Resiga et al., 2012) ont également présentéune méthodologie complète
pour l’évaluation et l’optimisation du comportement de l’écoulement à l’intérieur de l’aspirateur
d’une turbine hydraulique Francis, basée sur le modèle de vitesse d’entrée présenté dans (RF
Susan-Resiga, Muntean, Avellan, & Anton, 2011).

2.1. Paramétrage du profil de vitesse d’entrée de l’aspirateur

Dans ce projet, afin que la vitesse d'entrée puisse être décrite par une condition limite de profil de
vitesse, une hypothèse d'écoulement axisymétrique est posée pour l’écoulement à l'entrée de
l’aspirateur. Pour optimiser le profil de vitesse d’entrée de l’aspirateur, la première chose àfaire
est de trouver une façon de représenter le profil de vitesse mathématiquement. Essentiellement, les
profils de vitesse d’entrée de l’aspirateur sont des courbes qui décrivent les variations des
composantes de vitesse le long du rayon d’entrée de l’aspirateur. Comme pour les représentations
 10

numériques des courbes, il y a plusieurs méthodes permettant de paramétrer le profil de vitesse

d’entrée. Par exemple, la courbe peut être définie comme un ensemble de points d’interpolation en
utilisant différents types de splines ou définie par une fonction analytique, qui est obtenue en
ajustant une fonction àdes données expérimentales.

Figure 2.3 Système de coordonnées cylindriques pour la représentation des profils d’entrée

2.1.1 Modèle analytique

Le débit en aval d’une roue de turbine est un écoulement tourbillonnaire complexe. L’un des
modèles les plus simples et intuitifs pour la représentation de l’écoulement tourbillonnaire est le
modèle d’un corps solide en rotation, qui est la combinaison d’un profil de vitesse axiale uniforme
et un profil de vitesse angulaire uniforme. Mathématiquement, les profils de vitesse de rotation du
corps solide peuvent être définis comme suit :

Profil de vitesse axiale,

𝑉𝑎 (𝑟) = 𝑈0 2-3

Profil de vitesse tangentielle,

𝑉𝑡 (𝑟) = Ω0 𝑟 2-4

Où𝑉𝑎 est le profil de vitesse axiale, 𝑉𝑡 est le profil de vitesse tangentielle, 𝑟 est le rayon de la
section d'entrée, 𝑈0 est vitesse axiale et Ω0 est vitesse angulaire (voir Figure 2.3). Le profil de
vitesse radiale est déterminéàpartir du profil axial.
 11

Des études récentes ont montré que le profil de vitesse d’entrée de l’aspirateur peut être mieux
évalué et modélisé par une combinaison de tourbillons simples. Par l’analyse de données
expérimentales, Susan-Resiga et al. ont publiéun article dans lequel ils ont décrit que l’écoulement
en sortie d’une roue Francis peut être analytiquement représentépar l’ensemble suivant d’équations,
qui est une superposition de trois tourbillons élémentaires, àsavoir un corps solide en rotation, un
vortex de Batchelor en co-rotation et un vortex de Batchelor en contre-rotation (Romeo Susan-
Resiga et al., 2006).

Profil de vitesse axiale,

𝑟2 𝑟2
𝑉𝑎 (𝑟) = 𝑈0 + 𝑈1 exp (− ) + 𝑈2 exp (− ) 2-5
𝑅12 𝑅22

Profil de vitesse tangentielle,

𝑅12 𝑟2 𝑅22 𝑟2
𝑉𝑡 (𝑟) = Ω0 𝑟 + Ω1 [1 − exp (− 2 )] + Ω2 [1 − exp(− 2 )] 2-6
𝑟 𝑅1 𝑟 𝑅2

Où𝑉𝑎 est le profil de vitesse axiale, 𝑉𝑡 est le profil de vitesse tangentielle, 𝑟 est le rayon de la
section d'entrée, 𝑈0 , 𝑈1 , 𝑈2 sont les vitesses axiales, Ω0 , Ω1 , Ω2 sont les vitesses angulaires, 𝑅1 , 𝑅2
sont les rayons des cœurs des tourbillons.

Puis Gagnon et al. ont montréque les profils de vitesse axiale et tangentielle àla sortie d’une roue
Kaplan peuvent également être analytiquement approchés avec ces équations (J.-M. Gagnon et al.,
2012; J. Gagnon, Iliescu, Ciocan, & Deschênes, 2008).

Les modèles analytiques garantissent que les profils de vitesse optimaux trouvés par l’algorithme
d’optimisation sont physiques et sont plus susceptibles de pouvoir être générés par une roue de
turbine réelle. Cependant, de meilleures solutions pourraient exister et rester en dehors de l’espace
d’optimisation du modèle analytique.

2.1.2 Modèle d’interpolation

Une autre représentation qui pourrait être utilisée pour définir les profils de vitesse d’entrée de
l’aspirateur est le modèle d’interpolation, qui génère les profils de vitesse par interpolation d’un
ensemble de points de contrôle.
 12

Bien que toutes les techniques d’interpolation numériques puissent être utilisées comme un modèle
d’interpolation des profils de vitesse d’entrée, la spline cubique d’Hermite sera particulièrement
étudiée en raison de ses propriétés, qui incluent entre autres que :

- La courbe d’interpolation peut être contrôlée avec précision, car elle passe par tous les points de
contrôle.

- Les courbes sont interpolées par un polynôme cubique qui garantit la régularité de l’interpolation.

- Le polynôme cubique limite également l’influence des points de contrôle dans leurs intervalles
voisins.

Pour un ensemble donné de points 2D 𝑃⃗0 (𝑥0 , 𝑦0 ) … 𝑃⃗𝑛 et de vecteurs tangents correspondants
𝑃⃗′0 (𝑥 ′ 0 , 𝑦 ′ 0 ) … 𝑃⃗′𝑛 , l’interpolation entre 𝑃𝑘 et 𝑃𝑘+1 peut être donnée par

𝑃⃗𝑘
𝑃⃗
𝑃⃗𝑘 (𝑢) = [𝐻1 𝐻2 𝐻3 𝐻4 ] 𝑘+1
𝑃⃗′𝑘
[𝑃⃗ ′𝑘+1 ]

𝐻1 = 2𝑢3 − 3𝑢2 + 1 2-7

𝐻2 = −2𝑢3 + 3𝑢2
𝐻3 = 𝑢3 − 2𝑢2 + 𝑢
𝐻4 = 𝑢3 − 𝑢2

0≤𝑢≤1

L’interpolation monotone cubique d’Hermite est une variante des splines cubiques d’Hermite qui
préserve la monotonie de l’ensemble des données interpolées et évite les problèmes d’oscillation
de la spline cubique classique d’Hermite (Figure 2.4). De plus, elle fournit un moyen pour calculer
les vecteurs tangents aux points de contrôle, ce qui réduit de moitiéle nombre de variables de
conception de l’interpolation cubique d’Hermite classique et augmente considérablement
l’efficacité de l’optimisation.
 13

Figure 2.4 Comparaison entre la spline cubique monotone et classique d’Hermite

Dans l’interpolation monotone cubique d’Hermite (Fritsch & Carlson, 1980; Moler, 2008), pour
les points 2D, les vecteurs tangents 𝑃⃗ ′0 (𝑥 ′1 , 𝑦 ′1 ) … 𝑃⃗ ′𝑛 sont calculés avec les équations suivantes.

Calculez la pente des sécantes entre les points de contrôle voisins 𝑘 = 0, … , 𝑛 − 1

𝑦𝑘+1 − 𝑦𝑘
∆𝑘 =
𝑥𝑘+1 − 𝑥𝑘

Pour les points de contrôle internes 𝑘 = 1, … , 𝑛 − 1,

Si ∆𝑘−1 et ∆𝑘 ont le même signe,

∆𝑘−1 + ∆𝑘
𝑃⃗′𝑘 =
2

Si ∆𝑘−1 et ∆𝑘 ont des signes différents ou l’un d’entre eux est égal àzéro,

𝑃⃗′𝑘 = 0

Pour les points aux extrémités,

𝑃⃗′0 = ∆0

𝑃⃗ ′𝑛 = ∆𝑛−1
 14

2.2 Simulations CFD de l’écoulement interne dans un aspirateur

La simulation computationnelle en dynamique des fluides (Computational Fluid Dynamics ou
CFD) est une technique numérique pour simuler et analyser les problèmes de mécanique des fluides
en résolvant les équations de Navier-Stokes (N-S) àl’aide de méthodes numériques. Dans le cas
des simulations hydrauliques, le fluide est considérécomme incompressible et les équations de
Navier-Stokes dans cette forme sont,

∇∙𝑢
⃗ =0 2-8

𝜕𝑢
⃗ 1
+𝑢 ⃗ = − ∇p + 𝜈∇2 𝑢
⃗ ∙ ∇𝑢 ⃗ − ∇h 2-9
𝜕𝑡 ρ

Où𝑢
⃗ est le vecteur de vitesse, 𝜈 est la viscositécinétique et h est un terme source.

De nos jours, bien que les grappes de calcul àhaute performance soient beaucoup plus puissantes
que jamais, le coût de calcul de la résolution des équations de N-S dans des géométrie avec des
conditions limites et frontières complexes est encore énorme. Pour les problèmes d’optimisation
impliquant les simulations CFD, les défis deviennent encore plus sévères car la résolution d’un
problème d’optimisation nécessite d’évaluer les fonctions objectifs des centaines, voire des milliers
de fois, ce qui signifie l’exécution de la même quantité de simulations CFD. Pour éviter l’énorme
coût de calcul et pour rendre la méthodologie d’optimisation pratique, différents types d’équations
simplifiées ont étérésolus pour remplacer les équations complètes de N-S.

D’autre part, la simulation de l’écoulement à l’intérieur d’un aspirateur de turbine hydraulique est
un problème comportant de très grands défis pour les applications CFD puisque l’écoulement en
aval des roues est très fluctuant, tourbillonnant et non uniforme (Galván, Reggio, & Guibault,
2011). Un modèle trop simplifiéne peut pas reproduire correctement les structures turbulentes
complexes et le comportement de l’écoulement dans un aspirateur réel. Une méthode
d’optimisation, qui est construite en utilisant des simulations CFD trop simplifiées ne permet ainsi
pas d’obtenir des résultats utilisables.
 15

2.2.1 Modèles de basse fidélité

Les modèles mathématiques décrivant la physique de l'écoulement utilisés dans les simulations
CFD peuvent être classés en fonction de leur niveau de fidélité. Dans le cas de la simulation de
l’écoulement de l’aspirateur, qui a un grand nombre de Reynolds, l’écoulement est dominépar des
forces d’inertie et les effets visqueux sont significatifs seulement à l’intérieur des couches limites.
Par conséquent, les modèles à faible fidélité négligent l’effet visqueux du fluide, ce qui représente
la plus grande partie du coût de calcul de la simulation CFD, et ce qui permet, dans certains cas,
d’atteindre encore une précision acceptable. Dans les années 1980, lorsque la performance des
grappes d’ordinateurs n’était pas encore aussi avancée qu’aujourd’hui, ces modèles de basse
fidélité étaient de bons compromis entre la précision et l’efficacité. (Keck & Sick, 2008).
Aujourd’hui, avec l’utilisation d’ordinateurs modernes, les modèles àhaute-fidélitésont préférés
par les concepteurs de turbines, en raison de leur plus grande précision. Cependant, les modèles de
basse-fidélité ont toujours un avantage irremplaçable puisque le calcul peut être terminé en
quelques minutes seulement, voire même en quelques secondes. Ainsi, les modèles de basse-
fidélité sont fréquemment utilisés pour l’initialisation des modèles àhaute-fidélitépour améliorer
le taux de convergence de ces derniers. Dans le domaine de l’optimisation, les modèles de basse-
fidélitésont combinés avec le modèle àhaute-fidélitépour développer une méthode d’optimisation
multi-fidélité. Les modèles de basse-fidélité sont généralement utilisés comme fonction de
substitution pour filtrer les résultats prometteurs de sorte que le nombre de simulations des modèles
àhaute-fidélitépuisse être réduit.

2.2.1.1 Écoulement potentiel

Le modèle d’écoulement potentiel est l’un des premiers succès de la CFD dans la modélisation des
écoulements de turbines. Ce modèle suppose que l’écoulement est incompressible, non rotationnel
et non visqueux. Toutefois, il ne peut fournir un résultat raisonnablement précis pour la simulation
de la roue que lorsque l’état de fonctionnement est àproximitédu point de rendement optimal
(Keck & Sick, 2008).

L’écoulement potentiel est régi par les équations suivantes.

𝑢
⃗ = ∇𝜑 2-10
 16

∇2 𝜑 = 0 2-11

Où𝑢
⃗ est vecteur de vitesse et 𝜑 est le potentiel de vitesse.

Bahrami et al. ont proposé une méthodologie d’optimisation multi-fidélitéde conception des aubes
de roue de turbine, qui a utilisé un solveur d’écoulement potentiel comme substitut de basse-fidélité
au modèle àhaute-fidélité. Les études de cas ont montréque la méthode multi-fidélitéest capable
d’identifier une conception d’aube de turbine optimisée grâce àun effort de calcul relativement
faible(Salman Bahrami, Tribes, Devals, Vu, & Guibault, 2013; S Bahrami, Tribes, von Fellenberg,
Vu, & Guibault, 2014).

2.2.1.2 Équations d’Euler

Le modèle d’Euler est un autre modèle CFD de basse fidélité largement utilisédans les simulations
CFD àun stade précoce pour les turbines hydrauliques. Ce modèle est régi par les équations d’Euler
et suppose que l’écoulement est incompressible et non visqueux.

∇∙𝑢
⃗ =0 2-12

𝑑𝑢
⃗ 1
+𝑢
⃗ ∙ ∇𝑢
⃗ + ∇p = 0 2-13
𝑑𝑡 ρ

Où𝑢
⃗ est vecteur de vitesse et 𝑝 est la pression.

De la fin des années 1980 jusqu’au milieu des années 1990, les analyses basées sur les équations
d'Euler en 3D étaient l’un des outils CFD les plus populaires. Celles-ci ont fourni des résultats plus
précis que les modèles précédents, ce qui a permis aux chercheurs et aux concepteurs de turbines
de mieux comprendre le comportement de l’écoulement àintérieur des turbines Francis. (Keck &
Sick, 2008)

2.2.2 Modèles àhaute fidélité

Avec le développement rapide de la performance computationnelle, de plus en plus de modèles à

haute-fidélité ont été mis en œuvre pour la simulation des écoulements dans les turbines
hydrauliques. Ces modèles àhaute-fidélitéprennent en compte l’effet de la viscosité du fluide, ce
qui fournit des résultats beaucoup plus précis tout en augmentant significativement le coût de
 17

calcul. Les grappes actuelles d’ordinateurs à haute performance permettent aux chercheurs de
considérer les simulations basées sur les équations de Navier-Stokes en moyenne de Reynolds
(Reynolds-averaged Navier–Stokes ou RANS) en tant que candidat pour construire leurs boucles
d’optimisation à la place des modèles précédents de basse-fidélité.

Les équations RANS sont une formulation des équations de Navier-Stokes, moyennées sur une
certaine période de temps. Tous les champs physiques sont décomposés en deux parties : la partie
moyenne et la partie fluctuante. L’intégrale en temps de la partie fluctuante est supposée être égale
àzéro et seulement les valeurs moyennes du champ physique sont résolues.

∇∙𝑢
⃗ =0 2-14

𝑑𝑢
⃗ 1
+𝑢 ⃗ = − ∇p + 𝜈∇2 𝑢
⃗ ∙ ∇𝑢 ⃗ 2-15
𝑑𝑡 ρ

Où𝑢
⃗ est le vecteur de vitesse, 𝑝 est la pression et 𝜈 est la viscositécinématique.

Les équations RANS peuvent décrire les écoulements turbulents et réduire considérablement
l’énorme coût de calcul des équations complètes de Navier-Stokes. Cependant, le processus de
moyennage introduit des variables supplémentaires, les contraintes de Reynolds, qui rendent ce
système incomplet.

Pour fermer le système d’équations, on le complète en utilisant un modèle de turbulence, qui vient
avec de nouvelles hypothèses. Parmi les différents modèles de turbulence, les modèles à deux
équations sont très populaires pour la simulation des turbines hydrauliques. L’hypothèse
fondamentale qui permet de les développer suggère que les propriétés de l’écoulement, y compris
les contraintes de Reynold, sont isotropes. Cette hypothèse n'est plus valable dans le cas d’un
écoulement hautement turbulent, où le gradient de vitesse change comme par exemple dans un
aspirateur. Mais la résolution de ces équations reste le meilleur modèle qui puisse être utilisédans
une méthode d’optimisation où les simulations CFD sont utilisées comme fonctions objectif, en
raison des capacités computationnelles actuelles.
 18

2.2.2.1 Modèles de turbulence k-ϵ

Les modèles k-ε sont de loin les modèles de turbulence à deux équations les plus testés aujourd’hui
et largement utilisés dans l’industrie. Par conséquent, ils sont mis en œuvre dans la plupart des
codes CFD (Argyropoulos & Markatos, 2015).

Le modèle k-ε standard a été proposé d’abord dans (Launder & Spalding, 1974) et il suppose que
l’écoulement est entièrement turbulent et donc les effets de la viscositémoléculaire peuvent être
négligés. Deux nouvelles variables, l’énergie cinétique de turbulence k et son taux de dissipation ε,
et leurs équations de transport sont introduites pour fermer les équations RANS.

L’équation de l’énergie cinétique de turbulence (𝑘) :

𝜕𝑘 𝜕𝑘 𝜕 𝜇𝑡 𝜕𝑘 𝜕𝑢𝑖
𝜌 + 𝜌𝑢𝑖 = [(𝜇 + ) ] + 𝜏𝑖𝑗 − 𝜌𝜖 2-16
𝜕𝑡 𝜕𝑥𝑖 𝜕𝑥𝑗 𝜎𝑘 𝜕𝑥𝑗 𝜕𝑥𝑗

L’équation du taux de dissipation de l’énergie cinétique de turbulence (𝜖) :

𝜕𝜖 𝜕𝜖 𝜕 𝜇𝑡 𝜕𝜖 𝜖 𝜕𝑢𝑖 𝜖2
𝜌 + 𝜌𝑢𝑖 = [(𝜇 + ) ] + 𝐶1𝜖 𝜏𝑖𝑗 − 𝐶2𝜖 𝜌 2-17
𝜕𝑡 𝜕𝑥𝑖 𝜕𝑥𝑗 𝜎𝜖 𝜕𝑥𝑗 𝑘 𝜕𝑥𝑗 𝑘

Où 𝐶1𝜖 = 1,44 et 𝐶2𝜖 = 1,92 sont des constantes, 𝜎𝑘 = 1,0 , 𝜎𝜖 = 1,3 sont le nombre de
turbulence de Prandtl pour k et ϵ respectivement.

La viscositéturbulente 𝜇𝑡 est calculée par

𝑘2
𝜇𝑡 = 𝜌𝐶𝜇 2-18
𝜖

Avec 𝐶𝜇 = 0,09.

Grotjans et al. ont conclu que l’écoulement dans l’aspirateur était correctement simulé avec le
modèle k-ε (Grotjans, 2001).

Skåre et al. ont comparéles modèles k-ε standard et réalisable(Shih, Liou, Shabbir, Yang, & Zhu,
1994) pour la simulation de l’aspirateur de Kaplan et ont conclu que ces deux modèles ont prédit
des structures similaires et étaient seulement légèrement différents dans leur calcul des quantités
d’ingénierie (P Skåre, 2001).
 19

Petit et al. ont comparé les mises en œuvre du modèle k-ε standard entre CFX et OpenFOAM dans
la simulation de la turbine de l’aspirateur U9 Kaplan. Les résultats ont montréque les prédictions
de CFX et OpenFOAM sont presque identiques dans l’aspirateur (Petit, Nilsson, Vu, Manole, &
Leonsson, 2008)

Vu et al. ont proposéune méthodologie utilisant le modèle RANS et le modèle de turbulence k-ε
pour calculer l'écoulement en régime stationnaire et instationnaire dans un aspirateur coudéet ils
ont validéen comparant les résultats de simulation de CFX et OpenFOAM avec les résultats du
projet FLINDT (Vu, Devals, Zhang, Nennemann, & Guibault, 2011).

2.2.2.2 Modèles de turbulence k–ω

Les modèles k-ω sont également des modèles de turbulence largement utilisés dans le secteur
industriel. De nombreux scientifiques et ingénieurs ont contribué à ces modèles, mais le
développement le plus significatif a étéfait par (Wilcox, 1988) et la dernière version de ce modèle,
également connu sous le nom de modèle k-ω Wilcox (2006) (Wilcox, 2008), est présentécomme
suit :

L’équation de la viscositécinématique turbulente (𝜐𝑡 ) :

𝑘 2𝑆𝑖𝑗 𝑆𝑖𝑗 7
𝜐𝑡 = ̃ = max {𝜔, 𝐶𝑙𝑖𝑚 √
,𝜔 ∗
} , 𝐶𝑙𝑖𝑚 = . 2-19
𝜔
̃ 𝛽 8

L’équation de l’énergétique cinétique de turbulence (𝑘) :

𝜕𝑘 𝜕𝑘 𝜕𝑢̅𝑖 𝜕 𝑘 𝜕𝑘
+ 𝑢̅𝑗 = 𝜏𝑖𝑗 − 𝛽 ∗ 𝑘𝜔 + [(𝑣 + 𝜎 ∗ ) ] 2-20
𝜕𝑡 𝜕𝑥𝑗 𝜕𝑥𝑗 𝜕𝑥𝑗 𝜔 𝜕𝑥𝑗

L’équation du taux de dissipation spécifique (𝜔) :

𝜕𝜔 𝜕𝜔 𝜔 𝜕𝑢̅𝑖 𝜎𝑑 𝜕𝑘 𝜕𝜔 𝜕 𝑘 𝜕𝜔
+ 𝑢̅𝑗 = 𝛼 𝜏𝑖𝑗 − 𝛽𝜔2 + + [(𝑣 + 𝜎 ) ] 2-21
𝜕𝑡 𝜕𝑥𝑗 𝑘 𝜕𝑥𝑗 𝜔 𝜕𝑥𝑗 𝜕𝑥𝑗 𝜕𝑥𝑗 𝜔 𝜕𝑥𝑗

13 9 1 3 1 9
Avec 𝛼 = 25, 𝛽 = 𝛽0 𝑓𝛽 , 𝛽 ∗ = 100, 𝜎 = 2, 𝜎 ∗ = 5, 𝜎𝑑𝑜 = 8 , 𝛽0 = 125
 20

𝜕𝑘 𝜕𝜔
0, ≤0
𝜕𝑥𝑗 𝜕𝑥𝑗 1+85𝜒 Ω𝑖𝑗 Ω𝑗𝑘 𝑆̂𝑘𝑖 1 𝜕𝑢𝑚
𝜎𝑑 = { , 𝑓𝛽 = 1+100𝜒𝜔 , 𝜒𝜔 = | (𝛽∗𝜔) ̂
3 |, 𝑆𝑘𝑖 = 𝑆𝑘𝑖 − 2 𝜕𝑥 δki ,
𝜕𝑘 𝜕𝜔 𝜔 𝑚
𝜎𝑑𝑜 , >0
𝜕𝑥𝑗 𝜕𝑥𝑗

1 𝜕𝑢 𝜕𝑢 1 𝜕𝑢 𝜕𝑢 0 𝑖𝑓 𝑖 ≠ 𝑗
Ω𝑖𝑗 = 2 (𝜕𝑥 𝑖 − 𝜕𝑥𝑗), 𝑆𝑖𝑗 = 2 (𝜕𝑥𝑗 + 𝜕𝑥𝑗) et δij = {
𝑗 𝑖 𝑗 𝑖 1 𝑖𝑓 𝑖 = 𝑗

Pour de nombreux aspects, le modèle k-ω est mieux que le modèle k-ε standard. Il fournit une
prévision plus précise des couches limites soumises àdes gradients de pression négatifs. À cause
de cela, il est très appropriépour l’écoulement proche de la paroi, comme dans un aspirateur. Il ne
comprend également pas de fonctions d’amortissement non-linéaires qui sont nécessaires pour les
modèles k-ε.

Cependant, le modèle k-ω n’est pas aussi largement utilisés que les modèles k-ε, car il n’est pas
facile à converger. Dans de nombreux cas, le modèle k-ε standard est résolu comme une
initialisation pour améliorer la convergence de modèle k-ω, qui requièrent beaucoup de ressources
informatiques supplémentaires. Un autre problème bien connu du modèle k-ω est sa forte
sensibilité aux conditions limites de la quantité turbulente ω. Une différence significative de
prédictions de l’écoulement peut être observée si la valeur spécifiée de ω varie àl’entrée. Ce n’est
pas une bonne caractéristique de ce modèle de turbulence puisque généralement les conditions aux
limites de ω sont estimées plutôt que mesurées.

Liu et al. ont proposé une méthodologie d’optimisation pour la turbine Francis basée sur le modèle
de turbulence RNG k-ω (Liu, Wu, Wu, & Nishi, 2011).

2.2.2.3 Modèle k-ω de transport de contrainte de cisaillement (SST)

Le modèle k-ω de transport de contrainte de cisaillement (Shear Stress Transport Turbulence

Model ou SST) combine les modèles k-ε et k-ω classiques, qui sont deux modèles de turbulence
populaires. Le modèle k-ε est utilisé pour prédire l’écoulement libre et il est combinéau modèle k-
ω dans la partie proche de la paroi. Par conséquent, non seulement le problème bien connu du
modèle k-ω, àsavoir qu’il est trop sensible au niveau de turbulence de l’écoulement libre à l’entrée,
est-il évité, mais aussi le problème de l’inexactitude des modèles k-ε dans les couches limites est
également contourné. Le modèle SST donne des prédictions très précises de la séparation de
l’écoulement pour des gradients de pression adverses.
 21

Cependant, les modèles SST partagent les mêmes défauts que les modèles k-ω au point de vue de
la convergence.

L’équation d’énergie cinétique de turbulence (𝑘) :

𝜕𝑘 𝜕𝑘 𝜕 𝜕𝑘
+ 𝑢̅𝑗 = 𝑃𝑘 − 𝛽 ∗ 𝑘𝜔 + [(𝑣 + 𝜎𝑘 𝜈𝑇 ) ] 2-22
𝜕𝑡 𝜕𝑥𝑗 𝜕𝑥𝑗 𝜕𝑥𝑗

L’équation du taux de dissipation spécifique (𝜔) :

𝜕𝜔 𝜕𝜔 𝜕 𝜕𝜔 1 𝜕𝑘 𝜕𝜔
+ 𝑢̅𝑗 = 𝛼𝑆 2 − 𝛽𝜔2 + [(𝑣 + 𝜎𝜔 𝜈𝑇 ) ] + 2(1 − 𝐹1 )𝜎𝜔2 2-23
𝜕𝑡 𝜕𝑥𝑗 𝜕𝑥𝑗 𝜕𝑥𝑗 𝜔 𝜕𝑥𝑖 𝜕𝑥𝑖

Avec,

𝜙 = 𝜙1 𝐹1 + 𝜙2 (1 − 𝐹1 ) (𝜙 𝑒𝑠𝑡 𝛼, 𝛽, 𝜎𝑘 , 𝑒𝑡 𝜎𝜔 )
5 3 9
Où 𝛼1 = 9 , 𝛼2 = 0,44 , 𝛽1 = 40 , 𝛽2 = 0,0828 , 𝛽 ∗ = 100 , 𝜎𝑘1 = 0,85 , 𝜎𝑘2 = 1 , 𝜎𝜔1 = 0,5 , et

𝜎𝜔2 = 0,856.
4
√𝑘 500𝜈 4𝜎𝜔2 𝑘
𝐹1 = tanh {{min [max ( ∗ , 2 ), ]} }
𝛽 𝜔𝑦 𝑦 𝜔 𝐶𝐷𝑘𝜔 𝑦 2

𝜕𝑢𝑖
𝑃𝑘 = min (𝜏𝑖𝑗 ,10𝛽 ∗ 𝑘𝜔)
𝜕𝑥𝑗

1 𝜕𝑘 𝜕𝜔
𝐶𝐷𝑘𝜔 = max (2𝜌𝜎𝜔2 , 10−10 )
𝜔 𝜕𝑥𝑖 𝜕𝑥𝑖

Et la viscositécinématique turbulente est donnée par,

𝛼1 𝑘
𝜈𝑇 = 2-24
max(𝛼1 𝜔 , 𝑆𝐹2 )

2
2√𝑘 500𝜈
𝐹2 = tanh [[max ( ∗ , )] ]
𝛽 𝜔𝑦 𝑦 2 𝜔

H Nilsson et al. ont utiliséle modèle de turbulence Low-Re k-Omega SST pour mettre en œuvre
la simulation CFD pour l’aspirateur de la Turbine-99 (Nilsson & Cervantes, 2012)
 22

2.2.3 Traitement près de la paroi

La condition aux limites sur la paroi de l’aspirateur est une condition d'adhérence, ce qui implique
une vitesse nulle àla paroi, tandis que le fluide s’écoule àgrande vitesse ailleurs dans l’écoulement,
ce qui fait en sorte que la vitesse d’écoulement varie rapidement avec la distance normale dans la
zone proche de la paroi. Pour bien capturer le comportement de l’écoulement dans cette région, un
maillage très fin près de la paroi est nécessaire, ce qui implique une énorme quantitéde ressources
de calcul. Même en utilisant un maillage fin près de la paroi, il est bien connu que les modèles k-ε
ne peuvent fournir des résultats satisfaisants dans la région proche d’une paroi adhérante. Pour
éviter ces problèmes, les lois de paroi ont étéintroduites.

Les lois de paroi sont des modèles empiriques pour représenter l’écoulement proche de la paroi.
L’idée est de placer les premiers nœuds de calcul de la paroi à l’extérieur de la couche visqueuse.
Au lieu de simuler le comportement de l’écoulement à l’intérieur de la couche visqueuse en
résolvant le modèle de turbulence, les propriétés de l’écoulement àces nœuds sont directement
évaluées par un modèle empirique simple.

L’analyse expérimentale et dimensionnelle montre que le profil de vitesse moyen dans la zone
proche de la paroi est régi par la relation logarithmique suivante,

1
𝑢+ = ln(𝑦 + ) 2-25
𝜅

Avec,

𝑈𝑡
𝑢+ =
𝑢𝜏

𝜌∆𝑦𝑢𝜏
𝑦+ =
𝜇

𝜏𝜔
𝑢𝜏 = √
𝜌

𝑢+ est la vitesse adimensionnelle et 𝑦 + est la distance àla paroi adimensionnelle. 𝜏𝜔 est le

cisaillement àla paroi (wall shear stress) et 𝜅 est la constante de von Karman.
 23

Les lois de paroi sont généralement utilisées avec des modèles de turbulence k-ε pour éviter leur
principal inconvénient qui est l’inexactitude dans la couche limite. Elles sont également utilisées
avec d’autres modèles de turbulence, comme k-ω et SST, afin de réduire la densitédu maillage et
pour accélérer les simulations.
 24

2.3 Algorithmes d’optimisation numérique

Mathématiquement, pour un ensemble de paramètres de conception 𝑥, un problème d’optimisation
peut être défini comme,

min 𝐽(𝑥)

𝑡. 𝑞. 𝑔𝑗 (𝑥) ≤ 0 𝑗 = 1, … , 𝑚1
2-26
ℎ𝑘 (𝑥) = 0 𝑘 = 1, … , 𝑚2

𝑥𝑖 𝑙 ≤ 𝑥𝑖 ≤ 𝑥𝑖 𝑢 𝑖 = 1, … , 𝑛

𝐽(𝑥) est la fonction objectif, 𝑔𝑗 (𝑥) sont des contraintes d’inégalité et ℎ𝑘 (𝑥) sont des contraintes
d’égalité. On définit ainsi 𝑛 variables de conception 𝑥𝑖 avec leurs contraintes de bornes (𝑥𝑖 𝑙 and
𝑥𝑖 𝑢 ), 𝑚1 contraintes d’inégalité et 𝑚2 contraintes d’égalité.

Pour résoudre un problème d’optimisation continue avec des contraintes d’inégalité de ce genre, il
existe deux grandes catégories d’algorithmes d’optimisation, à savoir les méthodes basées sur les
gradients et les méthodes sans gradient.

2.3.1 Méthodes basées sur les gradients

Des algorithmes d’optimisation à base de gradient nécessitent l’existence de dérivées continues de
premier ordre de la fonction objectif et éventuellement des dérivées plus élevées. Par rapport aux
méthodes sans dérivées, elles ont généralement besoin d’un plus petit nombre de cycles de
conception pour les faire converger vers une solution optimale. Ainsi, l’algorithme basé sur les
gradients peut gérer àun grand nombre de variables d’optimisation.

Cependant, la convergence des méthodes basées sur les gradients est seulement garantie vers un
minimum local. Il est donc difficile d’obtenir des résultats satisfaisants dans un espace
d’optimisation bruité parce qu'il a beaucoup de minima locaux. De plus, ces méthodes ont
généralement de la difficulté lorsqu’elles traitent des fonctions objectifs qui sont discrètes,
discontinues, multimodales ou mélangées (discrètes et continues) (Nocedal & Wright, 2006).
 25

2.3.1.1 Calcul des dérivées

Les méthodes basées sur les gradients doivent calculer les dérivées de premier ordre de la fonction
objectif. D’une manière générale, cette information peut être évaluée par deux méthodes
différentes.

2.3.1.1.1 Analyse de sensibilitéanalytique

L’approche la plus puissante pour l’évaluation d’une dérivée de fonction est d’en déterminer son
équation analytique. Pour les fonctions objectif analytiques, l’équation analytique des dérivés peut
être obtenue avec un logiciel de calcul symbolique comme Mathematica ou MATLAB. Si la
fonction objectif est hautement non-linéaire et si son évaluation implique des méthodes itératives
complexes pour faire diminuer les résidus sous la tolérance spécifiée, comme dans le cas des
analyses structurelles et en simulation fluide, il est possible d’utiliser des logiciels de différentiation
automatique comme ADIC, ADIFOR et TAMC. Étant donné que toute fonction dans un
programme informatique peut être calculée par une série de séquences élémentaires impliquant la
multiplication, la division, l’addition, la soustraction et des fonctions trigonométriques et
exponentielles, ces logiciels peuvent analyser le code source du programme ligne par ligne, puis
produire la dérivée de toute opération par l’application de la règle de la dérivée en chaîne
automatiquement, sans aucune intervention de l’utilisateur ni la nécessité d’un graphe de
dépendance.

Les méthodes d’analyse sont très précises et efficaces. En général, l’évaluation des dérivés ne
nécessite que 4-5 fois supplémentaires le coût de calcul et l’occupation de la mémoire du code
original. Par conséquent, elles peuvent gérer un grand nombre de variables de conception, comme
par exemple pour des problèmes d’optimisation géométrique basée sur un ensemble de points de
contrôle.

Cependant, le code de source est obligatoire pour le processus de différenciation automatique.

Ainsi, ces méthodes ne peuvent pas s’appliquer aux codes commerciaux. Et la génération
automatique du code n’est pas parfaite ; pour une meilleure efficacité, du travail supplémentaire
doit être fait pour mettre en œuvre de meilleures structures de données et algorithmes.
 26

2.3.1.1.2 Approximation par différences finies

Si l’évaluation de la fonction objectif implique des codes commerciaux auxquels l’utilisateur ne

peut pas accéder ou des codes complexes àsources ouvertes qui sont très difficiles àmettre en
œuvre par la différenciation automatique, le seul choix devient alors l’approximation du gradient
par des méthodes numériques.

L’approximation par différences finies est une méthode universelle pour l’évaluation des dérivées
d’une fonction objectif, et elle est très facile à mettre en œuvre. Cependant, le coût de calcul de ce
procédéest proportionnel au nombre de variables de conception, et donc il ne permet pas de gérer
un problème d’optimisation avec un très grand nombre de variables de conception. De plus,
l’approximation par différences finies est difficile lors de la sélection d’un pas appropriée de
différenciation, en particulier pour un espace d’optimisation complexe et bruité.

2.3.1.2 Algorithme du gradient conjugué

Les algorithmes de gradients conjugués (GC) sont des algorithmes d’optimisation numérique à
base de gradient pour la résolution de problèmes d’optimisation non linéaires contraints (Nocedal
& Wright, 2006). Comme beaucoup de méthodes numériques, ils résolvent le problème
d’optimisation sans contrainte par une stratégie itérative qui commence par une estimation initiale
𝑥0 et qui ensuite met àjour le 𝑥 en effectuant des itérations de l’équation suivante :

𝑥𝑘+1 = 𝑥𝑘 + 𝛼𝑘 𝑝𝑘 2-27

Où𝑝𝑘 est une direction de recherche et 𝛼𝑘 est une taille de pas. Pour la première itération, la
direction de recherche donnée par l’inverse du gradient.

𝑝0 = −𝐽′ (𝑥0 ) 2-28

Pour les itérations suivantes, la direction de recherche 𝑝𝑘 peut être obtenue par l’algorithme de
Fletcher-Reeves, qui est un choix populaire pour le calcul de la direction de recherche parmi les
algorithmes de gradient conjugué.

𝐽′ (𝑥𝑘+1 )𝑇 𝐽′ (𝑥𝑘+1 )
𝑝𝑘 = 2-29
𝐽′ (𝑥𝑘 )𝑇 𝐽′ (𝑥𝑘 )
 27

La taille du pas 𝛼𝑘 peut être calculée par un algorithme de recherche de ligne (Nocedal & Wright,
2006).

Les méthodes de gradient conjuguén’exigent d’assemblée aucune matrice de calcul, réduisant ainsi
le stockage nécessaire, de sorte qu’elles ont de bonnes performances pour résoudre les problèmes
d’optimisation à grande échelle. Par contre, leur vitesse de convergence est plus lente que d’autres
algorithmes d’optimisation basés sur les gradients.

2.3.1.3 Algorithme de quasi-Newton

Les méthodes de quasi-Newton sont utilisées pour trouver les maxima et minima locaux des
fonctions d’optimisation (Nocedal & Wright, 2006). Il s’agit d’une alternative àla méthode de
Newton, qui est une méthode classique d’optimisation numérique itérative.

Dans la méthode de Newton, un point de départ 𝑥0 est donné. À toutes les itérations, une
approximation quadratique est construite via une expansion en série de Taylor vers la fonction
objectif àpartir du point 𝑥𝑘 actuel.

1
𝐽(𝑥) ≈ 𝐽(𝑥𝑘 ) + (𝑥 − 𝑥𝑘 )𝑇 𝐽′ (𝑥𝑘 ) + (𝑥 − 𝑥𝑘 )𝑇 𝐽" (𝑥𝑘 )(𝑥 − 𝑥𝑘 ) 2-30
2

La fonction d’approximation est minimisée pour obtenir le prochain point 𝑥𝑘+1 et la procédure est
répétée jusqu’à ce que le gradient de la fonction objectif converge vers un nombre suffisamment
petit.

𝑥𝑘+1 = 𝑥𝑘 − 𝐽" (𝑥𝑘 )−1 𝐽′ (𝑥𝑘 ) 2-31

La méthode de Newton a une très bonne performance pour la recherche locale et un taux de
convergence de 2ème ordre. Cependant, son coût de calcul peut être extrêmement coûteux si la
matrice hessienne de la fonction objectif 𝐽" (𝑥𝑘 ) provient d’une approximation par différences
finies.

Dans les méthodes de quasi-Newton, seules les dérivés de premier ordre 𝐽′ (𝑥𝑘 ) sont utilisées pour
se rapprocher de la matrice hessienne au lieu de la calculer àchaque itération. Les algorithmes
d’approximation de la matrice hessienne sont également disponibles pour les mises àjour de type
 28

quasi-Newton et l’algorithme de Broyden–Fletcher–Goldfarb–Shanno (BFGS) est le plus efficace

d’entre eux, de sorte qu’il est aussi le choix le plus populaire (Nocedal & Wright, 2006).

L’inverse de la matrice hessienne au point initial 𝐽" (𝑥0 )−1 pourrait être approchépar la méthode
des différences finies ou simplement approximée par une matrice identité. Puis les matrices
hessiennes suivantes sont obtenues par mise àjour par l’équation BFGS,

∆𝑔𝑘𝑇 𝐽" (𝑥𝑘 )−1 ∆𝑔𝑘 ∆𝑥𝑘 ∆𝑥𝑘𝑇

𝐽" (𝑥𝑘+1 )−1 = 𝐽" (𝑥𝑘 )−1 + (1 + ) 𝑇
∆𝑔𝑘𝑇 ∆𝑥𝑘 ∆𝑥𝑘 ∆𝑔𝑘
2-32
𝐽" (𝑥𝑘 )−1 ∆𝑔𝑘 ∆𝑥𝑘𝑇 + (𝐽" (𝑥𝑘 )−1 ∆𝑔𝑘 ∆𝑥𝑘𝑇 )𝑇
−
∆𝑔𝑘𝑇 ∆𝑥𝑘

Avec,

∆𝑔𝑘 = 𝐽′ (𝑥𝑘+1 ) − 𝐽′ (𝑥𝑘 )

∆𝑥𝑘 = 𝑥𝑘+1 − 𝑥𝑘

La méthode quasi-Newton avec mise àjour de BFGS réalise de très bonnes performances, même
pour les cas difficiles d’optimisation. Cependant, toutes les variantes de la méthode de Newton ou
quasi-Newton ont besoin de stocker et d’évaluer la matrice hessienne qui occupe beaucoup de
mémoire dans le cas de la résolution d’un système linéaire à grande échelle de sorte que
l’algorithme BFGS n’est guère utilisédans les problèmes d’optimisation à grande échelle.

2.3.2 Méthodes sans gradient

Les méthodes sans gradient peuvent trouver les solutions optimales en se basant entièrement sur
l’évaluation de la fonction objectif. Ainsi, les gradients et la matrice hessienne de la fonction
objectif n’ont pas àêtre calculés et stockés. Par rapport aux méthodes basées sur les gradients, les
méthodes sans gradient sont plus faciles à mettre en œuvre et plus universelles, et elles peuvent
résoudre différents types de problèmes d’optimisation.

Les méthodes sans gradient sont souvent utilisées dans les conditions suivantes (Conn, Scheinberg,
& Vicente, 2009),

- Les fonctions objectif non-différentiables ou bruitées et/ou contraintes, comprenant des minima
locaux multiples
 29

- Espace déconnecté, discret et/ou non-convexe

- Variables mixtes (discrètes, continues, permutation)

D’autre part, sans information sur les dérivés de la fonction objectif, le nombre de cycles de
conception des méthodes sans gradient est beaucoup plus grand que pour les méthodes basées sur
les gradients.

2.3.2.1 Algorithme génétique

Les algorithmes génétiques (AG) forment une grande famille d’approches d’optimisation qui ont
été inspirées par le processus de l’évolution naturelle des organismes, tout d’abord développépar
John Holland et ses collègues à la fin des années 1960 à l’Universitédu Michigan.

Les AG utilisent un codage binaire des variables de conception, aussi connu comme un
chromosome, au lieu des variables de conception elles-mêmes, ce qui leur permet de manipuler des
variables de conception continues et discrètes. Ils sont également capables de traiter des fonctions
objectif bruyantes ou non continues (Golberg, 1989).

Tous les AG sont basés sur trois opérateurs essentiels, àsavoir la sélection, le croisement, et la
mutation. L’opérateur de sélection est en charge du choix des chromosomes qui favorise la survie
du plus fort en donnant aux meilleurs chromosomes plus de chances de passer à la génération
suivante. L’opérateur de croisement met en œuvre un échange aléatoire de bits entre des paires
appariées de chromosomes, ce qui imite les processus de reproduction oùles traits génétiques sont
propagés. L’objectif principal de l’opérateur de croisement est de veiller à ce que la nouvelle
génération de chromosomes soit différente de la précédente. L’opérateur de mutation génère une
variation sur les chromosomes courants de façon aléatoire par retournement de bits, ce qui assure
la diversitédes chromosomes. Une procédure typique d’AG peut être décrite comme Figure 2.5

Contrairement aux méthodes basées sur les gradients, un AG commence ses recherches àpartir de
plusieurs points, et la population de chromosomes peut couvrir une large portion de l’espace
d’optimisation durant le processus d’exploration de sorte qu’il est un bon algorithme de recherche
global qui ne peut pas facilement être piégépar un minimum local tel que les méthodes basées sur
les gradients. De plus, la mise en œuvre de l’AG est simple et facilement parallélisable.
 30

Figure 2.5 Organigramme d’un algorithme génétique

Cependant, l’AG est une sorte de méthode à « force brute »et son taux de convergence est très lent
en comparaison avec d’autres algorithmes d’optimisation, en particulier avec les méthodes basées
sur les gradients. Et il a par ailleurs une faible capacitéde traitement des contraintes.

2.3.2.2 Recherche de patrons

La recherche de patrons est une famille d’algorithmes d’optimisation numérique sans gradients,
aussi connue comme méthode de recherche directe, qui calcule une séquence de points qui se
rapprochent d’un point optimal.

Figure 2.6 L'étape sondage de la recherche de coordonnées

 31

La recherche de coordonnées (Coordinate Search ou CS) est une variante rapide et simple dans la
famille de la recherche de patrons. Ces algorithmes commencent leurs recherches à partir d’un
point initial 𝑥0 . À chaque itération, l’algorithme modifie une variable de conception àla fois avec
la même taille de pas. Cette approche ressemble graphiquement à la recherche d’un ensemble de
points voisins autour d’un point actuel 𝑥𝑘 , avec une taille de pas fixe ∆𝑘 (Figure 2.6).

Si un meilleur point est trouvéparmi les points voisins, ce nouveau point devient le point courant
pour la prochaine itération. Lorsqu’il n’y a pas une telle augmentation ou diminution, la taille de
pas ∆𝑘 est réduite et la recherche est répété jusqu’à ce que les pas soient jugés suffisamment petits
(Figure 2.7).

Figure 2.7 Recherche de coordonnées

La recherche de patrons généralisée (Generalized Pattern Search ou GPS) est une généralisation de
la recherche de coordonnées (Figure 2.8). Par rapport aux méthodes de recherche de patrons
traditionnelles, une étape de recherche facultative est ajoutée au début de chaque itération. Dans
cette étape de recherche, un nombre fini de points de maillage peut être triéet évaluéde façon
opportuniste, ce qui améliore la capacité de recherche globale des méthodes traditionnelles de
recherche directe. Une stratégie populaire de recherche par étape est celle de la recherche le long
de la direction de sondage obtenue par l’étape de sondage précédente jusqu’à ce qu’un meilleur
point ne puisse être obtenu dans cette direction de sondage. Dans le cas de l’étape de mise àjour,
la taille du pas de la GPS est augmentée si un meilleur point a ététrouvélors des étapes précédentes
de cette itération. (Audet & Dennis Jr, 2002)
 32

Figure 2.8 Recherche de patrons généralisée

La recherche directe adaptative sur maillage (Mesh Adaptive Direct Search ou MADS) est une
généralisation de la GPS qui remplace le paramètre de taille de maillage ∆𝑘 de la GPS par un pas
𝑝
global ∆𝑚
𝑘 et ajoute un nouveau paramètre de taille de sondage ∆𝑘 . Deux tailles de maillages

différentes permettent àMADS de générer des directions de sondage plus souples que GPS (Figure
2.9).

Figure 2.9 Comparaison entre les directions de sondage de GPS (gauche) et MADS (droite)

Ainsi l’information des étapes de sondage peut être utilisée de manière plus adéquate. Il existe
plusieurs mises en œuvre de MADS qui définissent les différentes façons de générer des directions
de sondage telles que LT-MADS (Audet & Dennis Jr, 2006), QR-MADS (Van Dyke & Asaki,
2013) et OrthoMADS (Abramson, Audet, Dennis Jr, & Digabel, 2009).
 33

Wetter et al. ont comparéla recherche de patrons généralisée (GPS) et l’algorithme génétique (GA)
dans un problème d’optimisation de la conception d’un système de chauffage, ventilation et
climatisation (Heating, ventilation and air conditioning ou HVAC) comportant de grandes
discontinuités et plusieurs minima locaux dans la fonction objectif. Les résultats montrent que l’AG
a une capacitéde recherche plus globale que le GPS dans ce cas-ci (Wetter & Wright, 2003).

Bahrami et al. ont comparéles performances de la recherche directe sur maillage adaptatif (MADS)
et un algorithme évolutionnaire pour un problème d’optimisation d’une turbine Francis et ils ont
conclu que l’algorithme évolutionnaire a une capacité de recherche globale exceptionnelle tout en
sachant que MADS est plus performant dans la recherche locale. La capacitéde recherche globale
de MADS peut être améliorée en combinant l’échantillonnage hypercube latin pour obtenir les
points initiaux (Salman Bahrami, Tribes, von Fellenberg, Vu, & Guibault, 2015).

En conclusion, les algorithmes de recherche de patrons ont un meilleur taux de convergence et une
meilleure capacitéde recherche locale par rapport aux algorithmes génétiques. Cependant, ils ont
des performances relativement limitées pour ce qui est de leur parallélisme puisque les étapes de
recherche sont généralement des étapes effectuées séquentiellement et puisque seuls les étapes de
sondage peuvent être parallélisées. De plus, le nombre maximal de points voisins dans l’étape de
sondage dépend du nombre de variables de conception. Dans les algorithmes génétiques, le nombre
maximal d’évaluations simultanées est déterminé par la taille de la population de chromosomes.

2.4 Résumé
Ce projet vise à étendre les travaux de Galvan et propose de développer une méthodologie
d’optimisation pour déterminer les profils de vitesse d’entrée les plus appropriés pour les
aspirateurs àla sortie des roues de turbines Kaplan. Cette méthodologie sera mise en œuvre sur la
base de logiciels àcode source ouvert.

2.4.1 Paramétrage du profil de vitesse

La prioritédans ce projet est d’explorer des espaces d’optimisation aussi large que possible afin de
trouver des solutions optimales différentes de ce qu’aurait pu trouver un ingénieur en se basant sur
une approche traditionnelle. Un procédé de paramétrage de profil de vitesse flexible qui peut
représenter un espace d’optimisation aussi grand que possible est donc préférable. Ainsi, comparée
 34

àune représentation analytique, la représentation par interpolation des profils de vitesse offre une
plus grande flexibilité.

Un code python développéen interne par Andritz Hydro, appelédtpmaker, a étéchoisi comme
logiciel de paramétrage des profils de vitesse. Ce logiciel utilise la représentation par interpolation
pour contrôler le profil de vitesse et pour y intégrer un profil de couche limite àtravers une fonction
exponentielle près de la paroi. Le modèle de paramétrage du profil de vitesse dans dtpmaker permet
de créer un profil d’écoulement très flexible et d'introduire une couche limite physiquement réaliste
afin de combiner les avantages de la représentation par interpolation et de la représentation
analytique.

2.4.2 Simulation CFD

L’écoulement interne dans un aspirateur hydraulique est hautement turbulent et a une structure
complexe en raison de la non-uniformitéet des instabilités des champs de vitesse et de pression.
De ce fait, la prédiction de la performance de l’aspirateur avec une précision acceptable n’est pas
facile àobtenir sans l’aide de modèles de turbulence à haute-fidélitéet des maillages très fins. De
plus, l’algorithme d’optimisation permettra de tester un grand nombre de conditions aux limites
extrêmes sur le cas CFD. Alors le modèle CFD préférable doit offrir une bonne performance au
niveau de la convergence et de la robustesse.

D’autre part, les simulations CFD seront lancées des milliers de fois. Ainsi la différence de coût de
calcul entre les différents modèles CFD sera amplifiée et doit être prise en considération.

De l’examen précédent, deux types de modèles CFD peuvent être utilisés pour prédire l’écoulement
interne dans l’aspirateur. Dans le passé, les modèles de basse fidélitéétaient très populaires dans
le cas d’un problème d’optimisation avec la simulation CFD. Avec un ordinateur moderne, un cas
de simulation du modèle basse fidélitépeut être exécutéen quelques secondes. D’autre part, les
modèles de basse fidélité ignorent la simulation de l’effet de la viscosité et de la turbulence, ce qui
les rend imprécis dans leur évaluation des pertes d’énergie à l’intérieur de l’aspirateur.

Ces dernières années, le développement rapide des grappes de calcul à haute performance a
augmentéde manière significative les ressources de calcul disponibles, ce qui a fait en sorte que le
temps de calcul de la résolution des équations RANS chute àun niveau acceptable. Par conséquent,
les équations RANS sont choisies comme modèle CFD dans ce projet.
 35

En termes de modèle de turbulence, il y a quelques options qui se proposent. La revue précédente

montre que les modèles k-ε, les modèles k-ω et les modèles SST sont les choix les plus populaires
pour les simulations dans les turbines hydrauliques et qu’ils ont des caractéristiques différentes.

Du point de vue de la précision, les modèles k-ε fonctionnent bien pour les écoulements internes,
mais ne sont pas très précis lorsque que l’écoulement présente des gradients de pression
indésirables et une forte courbure à l’écoulement. Les modèles k-ω montrent une meilleure
précision que les k-ε modèles, en particulier dans les couches limites. Le modèle SST combine les
avantages des modèles k-ε et k-ω.

Du point de vue de la robustesse, les modèles k-ε sont les plus stables. Par rapport aux modèles k-
ε, les modèles k-ω sont plus sensibles à l’estimation initiale du champ d’écoulement. Les modèles
SST sont les plus sensibles. Ainsi les modèles k-ε ou k–ω ont généralement étérésolus en tant
qu’initialisation des modèles SST.

Vu le nombre de simulations, le modèle standard de turbulence k-ε a été choisi pour fermer des
équations RANS en raison de sa robustesse et sa bonne capacitéde convergence.

Une loi de paroi est également utilisée au mur, puisqu’elle réduit considérablement le nombre
d’éléments près des parois de l’aspirateur.

Pour profiter pleinement des avantages des grappes de calcul, des stratégies de parallélisation
doivent être introduites. Mais les codes CFD commerciaux coûtent énormément en frais de licence
lors de l’exécution de multiples cas de simulation en même temps, de sorte que leurs alternatives à
sources ouvertes deviennent un choix très attirant. Dans ce projet, OpenFOAM, un logiciel àcode
source ouvert très populaire, a étéchoisi comme le principal solveur CFD.

2.4.3 Algorithme d’optimisation

Dans ce projet, les variables de conception 𝑥𝑖 sont les variables de contrôle dans la méthodologie
de paramétrage du profil de la vitesse d’entrée et la seule fonction objectif est le facteur de perte
d’énergie dans l’aspirateur qui est évaluée par simulation CFD. La définition de la fonction objectif
est indiquée à l’équation 2-2 et le détail de son évaluation est donnée en 3.2.7. L’évaluation de la
fonction objectif implique le paramétrage du profil de vitesse et la simulation de l’écoulement. Par
conséquent, 𝐽(𝑥) est une fonction complexe et hautement non linéaire de 𝑥.
 36

Pour éviter des solutions peu pratiques, deux contraintes d’inégalité non linéaires devraient être
mises en œuvre. La première est une contrainte sur le nombre de tourbillon et la seconde une
contrainte sur le résidu final lors des simulations. La contrainte sur le nombre de tourbillon est
utilisée pour éliminer les solutions présentant un surplus de tourbillon puisque ce type
d’écoulement en aval indique une réduction de l’énergie extraite par la roue, entrainant une
diminution de l’efficacité globale de la turbine. La contrainte sur le résidu final vise àse débarrasser
des cas de simulation n’ayant pas bien convergé et dont la prédiction de la perte d’énergie dans
l’aspirateur n’est peut-être pas précise.

Pour résoudre un tel problème d’optimisation continue avec contraintes non-linéaires, plusieurs
algorithmes d’optimisation peuvent être choisis. Un algorithme d’optimisation approprié devra être
capable de faire face à des contraintes d’inégalité non linéaires et il doit avoir un bon taux de
convergence pour réduire les coûts et améliorer l’efficacité des processus itératifs de conception.

Les méthodes basées sur les gradients sont toujours le premier choix lorsque les dérivés de la
fonction objectif sont disponibles. Cependant, dans ce projet, plusieurs logiciels complexes
impliqués dans le processus d’évaluation de la fonction objectif impliquent un énorme travail pour
obtenir la solution analytique des dérivés. Si les dérivés étaient calculés avec des méthodes de
différences finies, le nombre d’évaluations exploserait, et atteindre un nombre comparable àcelui
des méthodes sans gradient. Et il est bien connu que la méthode basée sur les gradients pourrait
facilement rester piégée dans des minima locaux, en particulier dans le cas d’une fonction objectif
hautement non linéaire et non continue. Par conséquent, les méthodes sans gradient ont été
préférées dans ce projet.

Dans l’article de Galvan, l’algorithme d’optimisation est un algorithme génétique multi-île, qui est
une sorte d’algorithme évolutionnaire dont la force est la capacitéde recherche globale et le haut
degré de parallélisme. Cependant, le coût de calcul de l’algorithme évolutionnaire augmente
rapidement lorsque le nombre de variables de conception augmente et ainsi il n’est plus approprié
pour ce projet.

Dans ce projet, l’espace d’optimisation a été élargi en utilisant la représentation par interpolation
pour les profils de vitesse. Mais une autre hypothèse a étéposée, suggérant que les profils globaux
de vitesse optimale ne sont pas situés loin du profil de vitesse optimum local qui a ététrouvépar
la méthode traditionnelle utilisécomme le point de départ de l’optimisation dans ce projet. Un
 37

algorithme d’optimisation avec de bonnes performances de recherche locale a été préféré pour ce
type de problème d’optimisation.

Parmi les différents types de méthodes sans gradient, l’algorithme de recherche directe sur maillage
adaptatif (MADS), qui est un membre de la famille des algorithmes de recherche de patrons, a été
choisi pour sa bonne capacitéde recherche locale. Les détails de MADS ont étéexaminés àla
section 2.3.2.2.

NOMAD est une boîte à outils d’optimisation àcode source ouvert qui implémente l’algorithme
MADS pour résoudre des problèmes d’optimisation de type boîte noire sous des contraintes non
linéaires générales. NOMAD a étéchoisi comme logiciel d’optimisation dans ce projet en raison
de ses bonnes performances.
 38

CHAPITRE 3 MISE EN OEUVRE

Pour mettre en œuvre une méthodologie d’optimisation des profils de vitesse à l’entrée des
aspirateurs, un cadre logiciel appelécfdOpt a étédéveloppé. Python a étéchoisi comme langage
de programmation de cfdOpt en raison de sa diffusion et de sa popularitécroissante, qui s’explique
par l’intégration dans le langage de nombreux modules bien développés, rendant le développement
de nouvelles applications facile et rapide.

Figure 3.1 Diagramme de cfdOpt

Le cadre logiciel permet d’assembler une boucle d’optimisation comprenant trois parties
principales, à savoir le paramétrage, la simulation et l’optimisation. Au début de chaque boucle
d’optimisation, l’algorithme d’optimisation choisit un ou plusieurs ensembles de paramètres de
conception, en fonction de la description du problème d’optimisation àrésoudre. Les paramètres
de conception sont ensuite envoyés àla partie de paramétrage, dans laquelle ils sont convertis en
des conditions aux limites. La partie simulation prédit le facteur de perte d’énergie pour chaque
ensemble de conditions aux limites fourni par la partie de paramétrage en lançant des simulations
CFD. À la fin de la boucle d’optimisation, le facteur de perte d’énergie est retournéàla partie
d’optimisation en tant que valeur de la fonction objectif de l’optimisation. Ce processus est illustré
àla Figure 3.1.
 39

3.1 Paramétrage
Deux modèles de paramétrage ont été mis en œuvre dans cfdOpt, soit le modèle analytique et le
modèle par interpolation.

3.1.1 Modèle analytique

Pour fins de comparaison, la représentation analytique de l’écoulement en aval de la turbine de

Kaplan a été mise en œuvre en cfdOpt. Les détails de ce modèle analytique peuvent être trouvés à
la section 2.1.1.

3.1.2 Modèle d’interpolation

Nous détaillons ici la paramétrisation utilisant l’interpolation. Cette paramétrisation a étémise en

œuvre dans un programme Python, appelé“dtpmaker”, consituté d’un ensemble de scripts Python
développéà l’interne par Andritz Hydro et utilisépour générer des profils de vitesse réalistes pour
les paramètres de turbine donnés. Un profil typique de vitesse générépar dtpmaker comprend 3
segments différents, àsavoir le segment principal du profil, la couche aux limites intérieure et la
couche aux limites extérieure.

3.1.2.1 Principaux segments des profils

Les principaux segments des profils de vitesse sont les segments interpolés. Dans ce segment, les
profils de vitesse sont manipulés par un ensemble de points de contrôle. La mise en œuvre de la
méthode d’interpolation est basée sur SciPy.interpolate.PchipInterpolator, qui est une mise en
œuvre en Python des splines cubiques d’Hermite monotones. Les équations sont décrites à la
section 2.1.2.

Le programme dtpmaker permet aux utilisateurs de définir les coordonnées x et y pour les points
de contrôle utilisés pour l’interpolation. Dans ce projet, les coordonnées x des points de contrôle,
qui sont les distances radiales, sont fixées et elles sont réparties uniformément le long du rayon
d’entrée pour réduire le nombre de variables de conception.

Pour équilibrer le coût de calcul et la flexibilité dans la représentation des profils, 5 points de
contrôle sont utilisés pour représenter le profil de vitesse axiale et 4 pour le profil de vitesse
tangentielle (voir Figure 3.2). Le profil de vitesse radiale est donnépar l’équation suivante,
 40

𝑉𝑟 = 𝑉𝑎 ∗ sin(𝜃(𝑟)) 3-1

Où 𝜃(𝑟) est la combinaison linéaire de l’angle du moyeu et l’angle d’ouverture du cône de

l’aspirateur, lorsque l’on ne considère pas la présence du moyeu, l’angle au centre est nul.

3.1.2.2 Couches limites intérieures et extérieures

Pour les maillages de l’aspirateur sans moyeu àl’entrée, l’épaisseur du segment représentant la
couche limite intérieure a été mise à 0. La couche limite externe est contrôlée par un profil
représentépar une loi de puissance avec un exposant valant 1/7. L’épaisseur de la couche limite
extérieure peut être fixée pour chaque type de diffuseur.

Figure 3.2 Profil de vitesse typique pour les cas dtpmaker sans moyeu

3.1.2.3 Profil de vitesse axiale

Afin de comparer différents profils entre eux, le débit de l’écoulement d’entrée est fixé. Les points
de contrôle de la vitesse axiale peuvent être manipulés de façon àmettre à l’échelle le profil axial
pour atteindre un débit spécifié. Dans ce projet, la vitesse axiale moyenne a étéfixée àune valeur
typique de 3m/s pour tous les cas.
 41

La Figure 3.3 montre comment le segment principal du profil est corrigépour compenser la perte
d’écoulement dans la couche aux limites.

Figure 3.3 Profil de vitesse axiale corrigépour le débit

Figure 3.4 Profil de vitesse tangentielle (dtpmaker)

 42

3.1.2.4 Profil de vitesse tangentielle

Les points de contrôle de la vitesse tangentielle permettent de manipuler de façon précise ce profil
de vitesse. La Figure 3.4 illustre un profil typique et la plage de variation de points de contrôle
(flèches vertes). Le point de contrôle du profil de vitesse tangentielle au centre de l’aspirateur est
mis àzéro pour les cas sans moyeu central à l’entrée.

3.2 Simulation CFD

L'équation de Navier-Stokes pour un écoulement incompressible en moyenne de Reynolds a été
résolue avec la fermeture du modèle de turbulence k-ε standard et avec la méthode des volumes
finis. Deux codes CFD, OpenFOAM et CFX, ont été intégrés dans la méthodologie d’optimisation.
OpenFOAM était le principal solveur CFD et le support de CFX a été ajouté à des fins de
comparaison et de validation. Cependant, parce qu’OpenFOAM et CFX ont mis en œuvre leurs
solveurs différemment, plusieurs configurations doivent être effectuées afin de rendre les résultats
comparables. Dans OpenFOAM, la viscosité cinématique du fluide 𝜈 a été fixé à 8,92578 ∗
10−7 𝑚2 ∙ 𝑠 −1 pour correspondre àla valeur utilisée par CFX pour l’eau.

3.2.1 CFX

CFX est l’un des codes commerciaux de CFD parmi les plus populaires dans l’industrie et
largement utilisédans la simulation de turbomachines, de telle sorte que la simulation mise en place
et les résultats de CFX sont robustes et validés. Par conséquent, CFX a été choisi comme une
référence pour les calculs OpenFOAM dans ce projet.

3.2.2 OpenFOAM

OpenFOAM est un logiciel CFD àcode source ouvert qui a une grande base d’utilisateurs dans la
plupart des domaines de l’ingénierie et de la science. En tant que code CFD àcode source ouvert,
OpenFOAM est d’abord devenu populaire dans le domaine universitaire, car il permet aux
chercheurs de modifier et de contrôler tous les détails des simulations CFD. Après plusieurs années
de développement, OpenFOAM est devenu de plus en plus stable et dispose maintenant d’une vaste
gamme de caractéristiques qui attirent de plus en plus l’attention des utilisateurs industriels.
 43

OpenFOAM n’est pas aussi convivial pour les utilisateurs que ses homologues commerciaux et il
a une courbe d’apprentissage plus escarpée pour les débutants, car l’utilisateur doit comprendre et
ajuster tous les détails de ses simulations. Mais une fois qu’un cas standard a été mis en place, des
cas similaires peuvent être faits sans effort. Cette fonction a rendu OpenFOAM appropriépour les
chercheurs, qui se concentrent généralement sur des cas de simulation similaires, et les utilisateurs
d’optimisation, qui ont besoin d’exécuter des simulations similaires des milliers de fois pour
chercher la meilleure solution.

Il y a deux branches traditionnelles d’OpenFOAM, àsavoir la version de la Fondation OpenFOAM

et la version foam-extend. Dans ce projet-ci, la version foam-extend a étéchoisie pour sa vaste
gamme de caractéristiques et les bibliothèques de post-traitement qui sont utiles pour la mise en
place du cas d’optimisation.

3.2.3 Conditions aux limites

Les conditions aux limites pour la vitesse, la pression et la turbulence sont données ci-dessous.

Tableau 3.1 Conditions aux limites

Limites
Entrée Sortie Paroi
Champs

Vitesse (𝑢) [m/s] Profil 1D de vitesse Gradient nul Fixée à0

Fixer la valeur
Pression (𝑝) [Pa] Gradient nul Gradient nul
moyenne à0

Énergie cinétique de
Intensitéturbulente Gradient nul Gradient nul
turbulence (𝑘) [J/kg]

Taux de dissipation
de l’énergie cinétique Fonction de paroi
Longueur de mélange Gradient nul
de turbulence (𝜀) évolutive
[J/(kg·s)]

Les conditions aux limites d’entrée pour les vitesses sont définies par des profils 1D de vitesse, ce
qui est une condition aux limites axisymétrique et les valeurs aux nœuds de l’entrée sont calculées
 44

sur la base de la position radiale du nœud. Les conditions aux limites de sortie en vitesse sont
définies comme des gradients nuls, ce qui correspond àdes conditions aux limites de Neumann où
le gradient est fixé à 0. Cette condition aux limites suppose que le champ de vitesse soit
complètement développédans la section de sortie.

Pour la même raison, les conditions aux limites de pression sont définies comme un gradient nul à
l’entrée et àla paroi. À la sortie, la valeur moyenne de la pression a étéfixée àzéro Pa.

Dans le cas de l’énergie cinétique turbulente (𝑘), àla paroi et àla sortie, on impose un gradient nul
et àl’entrée une intensité turbulente égale à5% de la vitesse d'entrée, on calcule la condition aux
limites d’entrée pour 𝑘 par l’équation suivante.

3
𝑘= ∗ √𝑢𝑖𝑛𝑙𝑒𝑡 ∗ 𝐼𝑛𝑡𝑒𝑛𝑠𝑖𝑡é 𝑡𝑢𝑟𝑏𝑢𝑙𝑒𝑛𝑡𝑒 3-2
2

𝐼𝑛𝑡𝑒𝑛𝑠𝑖𝑡é 𝑡𝑢𝑟𝑏𝑢𝑙𝑒𝑛𝑡𝑒 = 5%

Pour la condition aux limites d’entrée du taux de dissipation de l’énergie cinétique turbulente (𝜀),
OpenFOAM et CFX ont mis en œuvre des modèles légèrement différents. Ainsi, des travaux
supplémentaires doivent être faits.

Pour CFX, la condition aux limites d’entrée du modèle de turbulence a été choisie en se basant sur
«l’échelle de l’intensité et de la longueur ». L’intensité de turbulence a été fixée à 5% et un
paramètre empirique d’échelle de longueur de Foucault a été utilisé. Le taux de dissipation de
l’énergie cinétique turbulente (𝜀) est donnépar les équations suivantes (Ansys, 2012),

𝑘1,5
𝜀= 3-3
É𝑐ℎ𝑒𝑙𝑙𝑒 𝑑𝑒 𝑙𝑜𝑛𝑔𝑢𝑒𝑢𝑟 𝑑𝑒 𝐹𝑜𝑢𝑐𝑎𝑢𝑙𝑡

É𝑐ℎ𝑒𝑙𝑙𝑒 𝑑𝑒 𝑙𝑜𝑛𝑔𝑢𝑒𝑢𝑟 𝑑𝑒 𝐹𝑜𝑢𝑐𝑎𝑢𝑙𝑡 = 0,005 [𝑚]

Dans OpenFOAM, les valeurs d’entrée pour 𝜀 sont données par un modèle de longueur de mélange.

𝐶𝜇0,75 ∗ 𝑘1,5
𝜀= 3-4
𝐿𝑜𝑛𝑔𝑢𝑒𝑢𝑟 𝑑𝑒 𝑚é𝑙𝑎𝑛𝑔𝑒

𝐶𝜇 = 0,09
 45

Pour correspondre àCFX, la longueur de mélange de OpenFOAM a étécalculée àpartir de

𝐿𝑜𝑛𝑔𝑢𝑒𝑢𝑟 𝑑𝑒 𝑚é𝑙𝑎𝑛𝑔𝑒 = 𝐶𝜇0,75 ∗ É𝑐ℎ𝑒𝑙𝑙𝑒 𝑑𝑒 𝑙𝑜𝑛𝑔𝑢𝑒𝑢𝑟 𝑑𝑒 𝐹𝑜𝑢𝑐𝑎𝑢𝑙𝑡 = 0,0008216

Les conditions aux limites du taux de dissipation de l’énergie cinétique turbulente (𝜀) ont également
étédéfinies par des gradient nuls sur la sortie et la paroi.

3.2.4 Traitement près de la paroi

Dans CFX, la paroi de l’aspirateur a étéfixée àla paroi adhérente, qui utilise les lois de paroi dans
le cas de tous les modèles de turbulence basés sur l’équation de ε (Ansys, 2012). Dans ce cas-ci, la
condition aux limites de l’énergie cinétique turbulente a étémise àzéro gradient et la loi de paroi
a calculé le taux de dissipation 𝜖, pour tous les premiers nœuds de la paroi de l’aspirateur en
utilisant les équations suivantes,

3⁄4
𝜌𝑢∗ 𝐶𝜇
𝜖= ∗ 𝑘 3⁄2 3-5
𝑦̃ 𝜅

𝑦̃ ∗ = max(𝑦 ∗ , 11,06)

𝜌𝑢∗ ∆𝑦
𝑦∗ =
𝜇
⁄
𝑢∗ = 𝐶𝜇1 4 𝑘 1⁄2

Où𝜌 est la densitédu fluide, 𝜇 est la viscositédu fluide, ∆𝑦 est la distance du premier nœud de la
paroi, 𝐶𝜇 est égal à0,09, 𝜅 est la constante de von Karman établie à0,41 dans CFX.

Dans OpenFOAM 2.3.0 (la version de la fondation) seule la fonction de paroi standard est mise en
œuvre. Cependant, dans la version foam-extend 3.2, une fonction de paroi évolutive, qui est
presque identique à celles de CFX, a été utilisée automatiquement avec le modèle standard de
turbulence k-ε. Seule une légère différence dans sa borne inférieure d’𝑦̃ ∗ .

𝑦̃ ∗ = max(𝑦 ∗ , 𝑦𝑙𝑏 )
 46

Dans foam-extend 3.2, 𝑦𝑙𝑏 est calculé par un processus itératif de l’équation suivante avec une
estimation initiale 𝑦𝑙𝑏 = 11

log(max(𝐸 ∗ 𝑦𝑙𝑏 , 1))

𝑦𝑙𝑏 =
𝜅

Dans le cas où𝜅 = 0,41 et 𝐸 = 9,8, 𝑦𝑙𝑏 vaut 11,53.

Bien que CFX15 et foam-extend 3.2 aient mis en œuvre la fonction de paroi évolutive
différemment, ils devraient fournir des résultats très proches.

3.2.5 Schéma de discrétisation

Dans le cas de CFX, le schéma àhaute résolution a étéutilisépour les termes d’advection et un
schéma d’ordre 1 a étéchoisi pour la turbulence numérique.

Pour assurer une meilleure précision, dans les cas d’OpenFOAM, le schéma de deuxième ordre a
étéchoisi.

3.2.6 Initialisation du champ d’écoulement

L’initialisation du champ d’écoulement affecte grandement la convergence des simulations CFD,
surtout quand la CFD est couplée avec un algorithme d’optimisation, ce qui génère un grand
nombre de conditions aux limites extrêmes pour la CFD.

Dans les cas de CFX, l’estimation initiale a été fixée à"Automatique" dans ce solveur. Dans cette
configuration, l’estimation initiale est donnée par l’interpolation linéaire entre les conditions aux
limites de l’entrée et de la sortie (Ansys, 2012). Dans le cas d’OpenFOAM, un écoulement potentiel
a étérésolu comme initialisation du modèle RANS.
 47

3.2.7 Évaluation de la fonction objectif

La fonction objectif est le facteur de perte d’énergie à l'intérieur de l’aspirateur,

1 1
𝑃 𝑑𝐴 −
𝐴𝑖𝑛 ∫𝑖𝑛 𝑡 𝐴𝑜𝑢𝑡 ∫𝑜𝑢𝑡 𝑃𝑡 𝑑𝐴
𝜁= 3-6
1 𝑄 2
2 𝐴𝑖𝑛 )
𝜌 (

1
𝐿𝑎 𝑚𝑜𝑦𝑒𝑛𝑛𝑒 𝑑𝑒 𝑙𝑎 𝑝𝑟𝑒𝑠𝑠𝑖𝑜𝑛 𝑡𝑜𝑡𝑎𝑙𝑒 à 𝑙 ′ 𝑒𝑛𝑡𝑟é𝑒: ∫ 𝑃 𝑑𝐴
𝐴𝑖𝑛 𝑖𝑛 𝑡

1
𝐿𝑎 𝑚𝑜𝑦𝑒𝑛𝑛𝑒 𝑑𝑒 𝑙𝑎 𝑝𝑟𝑒𝑠𝑠𝑖𝑜𝑛 𝑡𝑜𝑡𝑎𝑙𝑒 à 𝑙𝑎 𝑠𝑜𝑟𝑡𝑖𝑒: ∫ 𝑃𝑡 𝑑𝐴
𝐴𝑜𝑢𝑡 𝑜𝑢𝑡

1 𝑄 2
𝐸𝑛𝑒𝑟𝑔𝑖𝑒 𝑐𝑖𝑛é𝑡𝑖𝑞𝑢𝑒 𝑑′𝑒𝑛𝑡𝑟é𝑒 ∶ 𝜌 ( )
2 𝐴𝑖𝑛

Cette équation peut être facilement calculée àpartir des résultats de CFD en utilisant les équations

1 2 1 2
(𝑃𝑖𝑛 + 2 𝜌𝑈𝑖𝑛 ) − (𝑃𝑜𝑢𝑡 + 2 𝜌𝑈𝑜𝑢𝑡 )
𝜁= 3-7
1 2
2 𝜌𝑈𝑡ℎ𝑟𝑜𝑎𝑡

1
𝑃𝑖𝑛 = ∫ 𝑃𝑑𝐴
𝐴𝑖𝑛 𝑖𝑛

1
𝑃𝑜𝑢𝑡 = ∫ 𝑃𝑑𝐴
𝐴𝑜𝑢𝑡 𝑜𝑢𝑡

1
𝑈𝑖𝑛 = ∫ |𝑢⃗ | 𝑑𝐴
𝐴𝑖𝑛 𝑖𝑛

𝑄
𝑈𝑡ℎ𝑟𝑜𝑎𝑡 =
𝐴𝑖𝑛

𝑄
𝑈𝑜𝑢𝑡 =
𝐴𝑜𝑢𝑡

Ce facteur est très difficile àmesurer dans des essais réels avec des aspirateurs, àcause des limites
dans les techniques de mesure. Par conséquent, le facteur de perte d'énergie a été évalué
 48

approximativement dans les essais en utilisant des données provenant de points de mesure de la
vitesse et de la pression.

Pour cette raison, l’expression d’OpenFOAM pour l’évaluation du facteur de perte d’énergie a
également dû être modifiée pour mieux faire correspondre les résultats avec données de
l’expérience. (voir Annexe A)

L’évaluation du post traitement et de la fonction objectif a étéfaite à l’aide de swak4foam, une

bibliothèque au code source ouvert fournie avec OpenFOAM et étendant les fonctions de
manipulation d’expressions impliquant des objets d’OpenFOAM, et rendant plus facile pour les
utilisateurs la mise en œuvre des conditions initiales complexes, des conditions aux limites, et le
post-traitement.

La bibliothèque swak4foam permet d’évaluer la valeur de la pression totale à l’entrée et de

l’énergie cinétique d’entrée, directement, puisque les expressions sont évaluées sur le patch
d’entrée. Un patch est un ensemble de faces dans le maillage qui sont situées sur le contour du
maillage et leurs directions normales pointent toujours vers l’extérieur.

Figure 3.5 La partie d’extension du maillage de l’aspirateur

Toutefois, dans le cas de la pression totale de sortie, certains travaux supplémentaires doivent être
faits puisque les maillages de l’aspirateur ont une partie d’extension servant à éviter les
recirculations (voir Figure 3.5). Ainsi le patch de sortie n’est pas la position de sortie réelle de la
 49

géométrie de l’aspirateur. Pour évaluer la pression totale de sortie, OpenFOAM doit connaître les
faces exactes à l’intérieur du maillage de l’aspirateur pour évaluer l’expression de la pression totale
de sortie.

Il y a deux façons d’évaluer les expressions sur les faces à l’intérieur du maillage qui évaluent les
expressions sur une Surface définie ou sur une portion du maillage appelée faceZone.

3.2.7.1 Surface

La méthode la plus simple pour évaluer les expressions sur la sortie réelle de l’aspirateur est de
définir une surface sur la sortie en fournissant un point d’ancrage et un vecteur normal, puis
d’évaluer les expressions sur la surface ainsi définie. Cette approche est utile pour faire des
évaluations sur les géométries complexes ou les maillages non structurés, en raison de sa facilité
d’utilisation. Cependant, ce procédé introduit une erreur supplémentaire lors de la conversion entre
les données sur le maillage et la surface.

3.2.7.2 FaceZone

Dans OpenFOAM, un objet faceZone est un ensemble de faces dans le maillage comme un patch,
qui ne repose pas nécessairement sur les frontières du maillage. L’évaluation sur un faceZone est
considérée plus exacte que sur une surface puisque les expressions sont évaluées sur les faces
exactes des éléments qui sont sélectionnées par l’utilisateur.

Cependant, les utilisateurs doivent établir des critères pour sélectionner les faces spécifiques du
maillage. En général, les critères sont une gamme d’emplacements au centre de la face, ce qui
implique que les utilisateurs doivent connaître la position géométrique exacte des faces dans le
maillage. Pour des maillages structurés, la détermination des critères appropriés n’est pas difficile,
mais cela pourrait être très difficile dans le cas des maillages non structurés.

Les deux méthodes ont été mises en œuvre et comparées ; elles fournissent des résultats très proches.
Dans ce projet, la méthode faceZone a étéchoisie car les cas de simulation ont utilisédes maillages
structurés. De plus, dans certains cas, le facteur de perte d’énergie peut osciller périodiquement au
cours des itérations. Une méthode moyennage devrait être mise en œuvre pour une meilleure
précision de l’évaluation.
 50

3.2.8 Critères de convergence

Les simulations CFD visent à évaluer la fonction objectif pour l’optimisation, qui est la perte
d’énergie de l’écoulement à l’intérieur de l’aspirateur. Étant donnéque la fonction objectif est une
quantité moyenne du champ d’écoulement entier, la précision des détails de l’écoulement n’est pas
la principale préoccupation de la simulation. Les critères de convergence des simulations peuvent
donc être assouplis dans ce projet.

Figure 3.6 Le facteur de perte d’énergie (gauche) et les courbes résiduelles d’un cas typique de
CFD (droite)

Tel que montréàla Figure 3.6 pour un cas de simulation proche du point de meilleur rendement,
la courbe du facteur de perte d’énergie atteint l’état d’équilibre une fois que le résidu de la pression
est inférieure à10−3. Pour plus de sécurité, les critères de convergence de la simulation stipulent
que le résidu de la pression soit inférieur à10−5.

3.2.9 Validation de la cohérence entre OpenFOAM et CFX

3.2.9.1 Validation des codes de post-traitement

Les fonctions objectifs ont été évaluées avec la bibliothèque swak4foam dans les cas
d’OpenFOAM, et avec CFD-Post dans les cas CFX. En utilisant différents codes de post-
traitement, on pourrait introduire des erreurs supplémentaires. Ainsi, les résultats de la simulation
d’OpenFOAM dans le cas de l’aspirateur coudé (voir plus de détails à la section 4.2) ont été
exportés et post-traités dans CFD-Post à titre de comparaison. La comparaison a montré que
swak4foam et CFD-Post fournissent des résultats presque identiques (voir Figure 3.7) ; la majorité
de la différence entre OpenFOAM et CFX vient de leurs solveurs CFD.
 51

Il vaut la peine de mentionner qu’il y a différentes façons d’exporter les résultats de simulations
d’OpenFOAM vers CFD-Post pour le post-traitement. OpenFOAM fournit deux utilitaires appelés
foamToFluentData et foamToCGNS qui permettent aux utilisateurs de convertir leurs solutions
d’écoulement au format correspondant. Cependant, foamToFluentData ne peut pas convertir
correctement la distribution de vitesse de la condition aux limites spécifiée comme un
profile1DfixedValue dans OpenFOAM, qui est la condition aux limites d’entrée utilisée avec
OpenFOAM dans ce projet. À l’heure actuelle, foamToCGNS est la seule façon de faire le travail.

Figure 3.7 Comparaison entre swak4foam et CFD-Post

3.2.9.2 La validation des conditions aux limites d’entrée

Étant donné que les conditions aux limites d’entrée des modèles de turbulence sont différentes dans
les cas d’OpenFOAM et CFX, les profils d’entrée de 𝑘 et 𝜀 utilisés dans les deux logiciels ont été
visualisés et comparés pour valider leur cohérence.

À titre de comparaison, un profil de vitesse d’entrée typique (Figure 3.8) a étéchargéàla fois dans
OpenFOAM et CFX. Puis les profils d’entrée de 𝑘 et 𝜀 ont étéévalués par CFD-Post qui est l’outil
de post-traitement pour CFX. Les résultats montrent que les différentes conditions aux limites
d’entrée en OpenFOAM et CFX fournissaient des conditions aux limites similaires mais non
parfaitement identiques pour le modèle de turbulence 𝑘 − 𝜀 (Figure 3.9 et Figure 3.10).
 52

Figure 3.8 Le profil de vitesse d’entrée pour la comparaison des conditions aux limites d’entrée

Figure 3.9 Comparaison des profils d’entrée 𝑘

 53

Figure 3.10 Comparaison des profils d’entrée 𝜀

3.3 Optimisation

3.3.1 Point de départ de l’optimisation

L’algorithme MADS a une très bonne capacité de recherche locale et, dans ce projet, on a supposé
que la solution optimale n’était pas située loin du profil de vitesse optimal local qui a étédéterminé
par le test de rotation du corps solide (Solid Body Rotation ou SBR).

L’écoulement SBR est un écoulement tourbillonnaire avec une vitesse axiale constante et une
vitesse angulaire constante. La Figure 4.10 est un exemple typique de profils de vitesse SBR avec
ajout de couche limite. Dans ce projet, tous les maillages ont étégénérés par un générateur de
maillage structuré, spécifiquement développépour mailler les aspirateur appeléDTMesh. Tous les
diamètres d’entrée ont été mis à l’échelle de 1m. Par conséquent, un test SBR standard dont le
profil de vitesse axiale est de 3m/s et un faible nombre de tourbillon variable de 0 à0,5 a étéutilisés.
La Figure 4.3 illustre un nombre de tourbillon typique – la courbe du facteur de perte d’énergie est
obtenue àpartir d’un balayage en nombre de tourbillon pour des profils SBR. Pour l’écoulement
SBR, le nombre de tourbillon peut être évaluéde façon approximative par l’équation simplifiée
suivante.
 54

max(𝑉𝑡 )
𝑆= 3-8
2Va

De la courbe illustrée àla Figure 3.7, une conclusion qui peut être obtenue est que l’efficacité de
l’aspirateur diminue si l’écoulements d’entrée n’a pas assez ou trop de tourbillon. Il existe un point
de rendement optimal qui permet à l’aspirateur d’extraire le plus d’énergie à partir de l’écoulement
courant àla sortie de la turbine.

3.3.2 Contraintes non-linéaires

3.3.2.1 Contrainte sur le nombre de tourbillon

Par définition, le nombre de tourbillon est le rapport de la quantitéde mouvement de rotation de

l’écoulement axial sur la quantité de mouvement de rotation de l’écoulement dans la direction
tangentielle. Du point de vue de la conception de la turbine, les concepteurs préfèrent que
l’écoulement en aval de la turbine soit non-tourbillonnant, ce qui signifie que la majeure partie de
l’énergie cinétique de tourbillon dans l’écoulement d’entrée ait été extraite et transformée en
énergie cinétique de rotation de la turbine.

Par conséquent, il est nécessaire de mettre en œuvre une contrainte d’optimisation qui tente
d’éliminer les conditions d’entrée qui génèrent un écoulement ayant un très grand nombre de
tourbillon à l’entrée de l’aspirateur. En se basant sur l’analyse du comportement de plusieurs
aspirateurs, une valeur maximale de 0.4 pour le nombre de tourbillon a été choisie pour cette
recherche, et une contrainte a été ajoutée dans la définition du problème d’optimisation afin
d’éliminer immédiatement les solutions trop éloignées de la solution optimale recherchée. La
contrainte sur le nombre de tourbillon s’exprime ainsi comme :

𝑆 ≤ 0.4 3-9

3.3.2.2 Contrainte sur le résidu final

Certaines conditions de profils d’entrée extrêmes et irréalistes peuvent provoquer un écoulement

instable ou périodique à l’intérieur de l’aspirateur. Ces conditions doivent ensuite être évalués à
l’aide de simulation CFD au cours du processus d’optimisation. Ces écoulements instables peuvent
atteindre de très petits facteurs de perte d’énergie dans les évaluations, mais ces résultats doivent
 55

être considérés comme inexacts puisque les simulations en régime permanent ne sont pas
complètement convergées. Si ces résultats erronés sont retournés à l’algorithme d’optimisation, la
solution optimale souhaitée peut ne pas être atteinte. Par conséquent, une contrainte sur les résidus
finaux des simulations CFD a été mise en œuvre pour les éliminer.

Figure 3.11 Courbe résiduelle typique d’OpenFOAM

D’après les résultats des simulations avec OpenFOAM, la valeur du résidu final du champ de
pression a étéchoisie comme critère de la contrainte sur le résidu final puisque le champ de pression
possède le plus haut résidu parmi toute les variables dans la plupart des cas d’optimisation. De plus,
une étape de moyennage a été introduite afin de réduire l’impact de l’oscillation des résidus (Figure
3.11). Typiquement, une moyenne de 300 itérations était utilisée. Dans les cas d’optimisation réels,
les cas d’OpenFOAM dont le résidu moyen de pression finale est plus élevéque 1e-4, ont été
ignorés.

Résidu final ≤ 10−4 3-10

3.3.3 NOMAD et parallélisation

Comme la plupart des codes de calcul scientifique, OpenFOAM et CFX fonctionnent en parallèle
sur des processeurs distribués avec l’interface standard de passage de messages (Message Passing
Interface ou MPI). Cependant, en raison des défauts des méthodes multi-grilles (Chow, Falgout,
Hu, Tuminaro, & Yang, 2006), qui sont utilisées pour résoudre des systèmes linéaires en CFD,
 56

l’efficacité parallèle globale de certains cas CFD est limitée. Par conséquent, la mise en parallèle
de l’algorithme d’optimisation est également importante.

NOMAD est une mise en œuvre C ++ de l’algorithme Mesh Adaptive Direct Search (MADS) de
Audet et Dennis (Audet & Dennis Jr, 2006), qui est une méthode de recherche de patrons sur la
base du calcul non lisse de Clarke (Clarke, 1990). Il peut être utilisépour résoudre des problèmes
d’optimisation de type boîte noire avec des contraintes non linéaire. Une version MATLAB
développée par Mark Abramson est également disponible.

Par défaut, NOMAD fonctionne selon un mode d’exécution en série (illustréàla Figure 3.12), qui
a une meilleure compatibilitéavec les logiciels d’évaluation de la fonction objectif et est plus facile
à mettre en œuvre. Ce mode est adapté pour les problèmes d’optimisation à petite échelle ou pour
des exécutions sur un ordinateur personnel.

Figure 3.12 Mode série NOMAD

Cependant, les scénarios cfdOpt ont été testés sur Briaree, une grappe à haute performance de
Calcul Québec. Pour en tirer profit, le parallélisme au niveau de l’algorithme d’optimisation
devient très important.
 57

3.3.3.1 Mode MPI

NOMAD supporte également l’évaluation en parallèle des fonctions objectif en utilisant MPI. En
mode MPI, les processus multiples de NOMAD fonctionnent simultanément et chaque processus
gère une simulation CFD individuelle de sorte que tous les cas CFD peuvent être lancés et terminés
séparément. Cette caractéristique permet de mettre en œuvre une stratégie d’évaluation par
NOMAD plus souple, comme la recherche opportuniste, qui commence une nouvelle itération
d’optimisation aussitôt qu’une nouvelle meilleure solution a été trouvée, plutôt que d’attendre
jusqu’à ce que toutes les simulations CFD terminent.

Cependant, le mode MPI de NOMAD engendre des problèmes de compatibilité lorsque les
évaluations de la fonction objectif impliquent OpenFOAM, dont la parallélisation est également
basée sur le cadre MPI (illustréàla Figure 3.13).

Figure 3.13 Mode MPI NOMAD

3.3.3.2 Mode d’évaluation en blocs

En mode d’évaluation en blocs, NOMAD fonctionne en série, mais il fournit plus d’un ensemble
de variables de conception d’optimisation en une seule fois et il reçoit plus d’un ensemble de
résultats d’évaluations de la fonction objectif (illustréàla Figure 3.14). Les utilisateurs doivent
donc développer un nouveau niveau appelé ‘wrapper’(emballeur) d’optimisation pour gérer la
 58

répartition des tâches de simulation et la collecte des résultats de simulation. Puisque NOMAD
fonctionne en mode série, il n’y a plus de problèmes de compatibilité avec d’autres programmes
parallèles à base de MPI, comme OpenFOAM. D’autre part, dans ce mode, si NOMAD fournit
plus d’un ensemble de variables de conception, il ne peut pas recevoir les résultats de la simulation
avant que toutes les évaluations des variables de conception soient effectuées.

Figure 3.14 Mode d’évaluation de bloc NOMAD

Pour éviter les problèmes de compatibilité, NOMAD a étéréglé sur le mode d’évaluation en blocs
dans ce projet.
 59

CHAPITRE 4 RÉSULTATS

4.1 Cas du diffuseur conique

Le diffuseur conique constitue, au niveau géométrique, le modèle d’aspirateur le plus simple que
l’on peut concevoir puisqu’il est décrit àl’aide de seulement 3 paramètres de conception, àsavoir
le rayon d’entrée, l’angle d’ouverture et le taux d’expansion de la zone de diffusion. Par
conséquent, cette géométrie est très appropriée pour la validation et la comparaison de différents
modèles CFD. Dans ce projet, un cas de diffuseur conique est utilisépour calibrer la configuration
d’OpenFOAM sur celle de CFX, afin d’obtenir des résultats précis de simulation.

La géométrie du diffuseur conique utilisée est représentée àla Figure 4.1. Le diamètre d’entrée du
diffuseur conique est normaliséà1m, et le ratio d’aire entre l’entrée et la sortie de la zone est fixé
à 4. Une partie d’extension a été ajoutée àla sortie de l’aspirateur pour éviter les recirculations,
dont la longueur est de 3 fois le diamètre d’entrée.

Figure 4.1 Dimensions du diffuseur conique

 60

4.1.1 Maillage

Les maillages du diffuseur conique utilisés dans ce projet ont étégénérés par un code de génération
de maillage interne appeléDTMesh. Trois maillages structurés différents ont étégénérépour la
même géométrie de diffuseur conique pour vérifier la dépendance au maillage.

Figure 4.2 Maillage de diffuseur conique (Moyen)

Tableau 4.1 Information concernant les maillages du diffuseur conique

Densitéde maillage Grossier Moyen Fin

Nombre de nœuds 58536 272023 450296

Pour comparer les résultats numériques, les courbes de facteur de perte d’énergie en fonction du
nombre de tourbillon ont ététracées pour les trois maillages avec OpenFOAM (Figure 4.3) et CFX
(Figure 4.4) respectivement.

Dans le cas d’OpenFOAM, les maillages fins et moyens fournissent des résultats très similaires
pour la prévision du rendement àcharge partielle, tandis que les résultats sur le maillage grossier
sont plus éloignés. Les trois maillages prédisent des valeurs très proches pour l’efficacitéàpleine
charge.
 61

Figure 4.3 Courbes de facteur de perte d’énergie en fonction du nombre de tourbillon

(OpenFOAM)

Figure 4.4 Courbes de facteur de perte d’énergie en fonction du nombre de tourbillon (CFX)
 62

Dans le cas de CFX, la prédiction de l’efficacité globale avec les maillages moyen et fin est très
semblable. Cependant, les résultats obtenus avec le maillage grossier s’éloignent des résultats
obtenus sur les maillages plus raffinés.

Figure 4.5 Courbes de facteur de perte d’énergie et nombre de tourbillon (OpenFOAM et CFX)

La comparaison entre les prévisions de CFX et d’OpenFOAM (Figure 4.5) montrent un bon accord
pour des conditions proche du point de meilleur rendement (BEP). Cependant, la différence
augmente pour les conditions de charge partielle et de pleine charge. Ces écarts peuvent être
causées par plusieurs facteurs, àsavoir une mise en œuvre différente du modèle de turbulence, un
choix différent dans les méthodes de résolution et certaines difficultés de convergence pour des
conditions loin du BEP.

Il est également nécessaire de vérifier les erreurs introduites en utilisant différents outils de post-
traitement étant donné que les facteurs de perte d’énergie des simulations d’OpenFOAM et de CFX
ont étéévalués par swak4foam et CFD-Post, respectivement. Toutes les solutions d’OpenFOAM
ont été exportées vers CFD-Post pour un second post-traitement, et la comparaison entre les
résultats de CFD-Post et swak4foam est présentée àla Figure 4.6. La comparaison montre que ces
deux outils fournissent des résultats très similaires et la différence entre les cas OpenFOAM et
CFX proviennent davantage de différences entre les solveurs plutôt que des outils de post-
traitement.
 63

Figure 4.6 Comparaison entre les différents outils de post-traitement

Les courbes de convergence du résidu ont ététracées pour tous les cas de simulation afin d’observer
leur régularité. Dans le cas d’OpenFOAM, les cas de simulation dont le nombre de tourbillon est
inférieur au BEP ne peuvent pas atteindre les critères de convergence avant le dernier intervalle de
temps (Figure 4.7) tandis que tous les cas de simulation dont le nombre de tourbillon est plus élevé
que le BEP ont atteint les critères de convergence et se sont donc arrêtées avant la limite sur le
nombre d’itérations (Figure 4.8).

Dans le cas de CFX, la faible densitédu maillage a une meilleure chance de converger avant le
dernier intervalle de temps. Les informations détaillées de la convergence sont présentes dans le
Tableau 4.2.
 64

Figure 4.7 Courbes résiduelles typiques pour les cas de non-convergence dans CFX (gauche) et
OpenFOAM (droite)

Figure 4.8 Courbes résiduelles typiques pour les cas de convergence CFX (gauche) et
OpenFOAM (droite)
 65

Tableau 4.2 Convergence de tous les cas de simulation dans la Figure 4.5

(✗ - les cas de non-convergence, ✓ - les cas de convergence)

Solveur OpenFOAM CFX

Maillage

Tourbillon Grossier Moyen Fin Grossier Moyen Fin

0 ✗ ✗ ✗ ✗ ✓ ✗

0.0404 ✗ ✗ ✗ ✗ ✗ ✗

0.0807 ✗ ✗ ✗ ✗ ✗ ✗

0.1211 ✗ ✗ ✗ ✗ ✗ ✗

0.1615 ✗ ✗ ✗ ✗ ✗ ✗

0.2019 ✗ ✗ ✗ ✗ ✗ ✗

0.2422 ✗ ✗ ✗ ✗ ✗ ✗

0.2826 ✗ ✗ ✗ ✗ ✗ ✗

0.3230 ✗ ✗ ✗ ✓ ✗ ✗

0.3633 ✗ ✗ ✗ ✓ ✗ ✗

0.4037 ✓ ✓ ✓ ✓ ✗ ✗

0.4441 ✓ ✓ ✓ ✓ ✓ ✓

0.4844 ✓ ✓ ✓ ✓ ✓ ✓
 66

Dans les cas CFX, le maillage grossier a une meilleure convergence que les maillages moyen et
fin. Cependant, concernant l’aspect de la prédiction du coefficient de perte d’énergie, les maillages
moyens et fins donnent des résultats presque identiques. Dans les cas d’OpenFOAM, les maillages
moyens et fins donnent également des résultats très proche. Par conséquent, le maillage moyen a
été choisi pour être utilisé dans l’optimisation.

4.1.2 Solution initiale pour l’optimisation

Les profils de vitesse d’un corps solide en rotation au BEP prédits par OpenFOAM , ont étéchoisis
comme point de départ pour l’optimisation (Figure 4.9). Les profils de vitesse au BEP ont été
approchés par dtpmaker et tracés àla Figure 4.10. La valeur du facteur de perte d’énergie pour ce
point d’opération est 0,3534.

Figure 4.9 Courbe du facteur de perte d’énergie en fonction du nombre de tourbillon (Diffuseur
conique)
 67

Figure 4.10 Point de départ de l’optimisation

4.1.3 Résultats de l’optimisation

Le cas de l’optimisation du diffuseur conique vise à vérifier l’exactitude de la mise en œuvre de

cfdOpt et a montré la cohérence entre les résultats d’OpenFOAM et de CFX (Figure 4.11). Dans
ce cas-ci, l’épaisseur de la couche limite extérieure a étéfixée à5% du diamètre dans dtpmaker,
selon les pratiques recommandées par Andritz pour le maillage du diffuseur conique. L’angle
d’ouverture a été fixéà12 degrés et tous les paramètres de dtpmaker ont gardéleurs valeurs de
défaut. La plage de variation de variables de conception a étéfixée entre 0 et 1.
 68

Figure 4.11 Le graphique de convergence NOMAD (Diffuseur conique)

Figure 4.12 Solution optimale du diffuseur conique

La solution optimale du cas du diffuseur conique est représentée àla Figure 4.12. Le facteur de
perte d’énergie de la solution est de 0,036798. Par rapport au point de départ, dont le facteur de la
perte d’énergie est 0,3534, la solution optimale réduit de près de 90% la perte d’énergie à l’intérieur
 69

du diffuseur. Les Figure 4.13 et Figure 4.14 montrent la convergence de la simulation CFD de la
solution optimale avec OpenFOAM.

Figure 4.13 Courbe de facteur de perte d’énergie (Diffuseur conique)

Figure 4.14 Courbes résiduelles (Diffuseur conique)

 70

Le résultat de la simulation du point optimal a étévisualisécomme Figure 4.15 et Figure 4.16.

Figure 4.15 Champ de vitesse de la solution optimale (OpenFOAM)

Figure 4.16 Contour de la pression de la solution optimale (OpenFOAM)

 71

Le résultat de la simulation est vérifiépar CFX et visualisécomme Figure 4.17 et Figure 4.18. Le
facteur de perte d’énergie prédit par CFX est 0,0397682.

Figure 4.17 Champ de vitesse de la solution optimale (CFX)

Figure 4.18 Contour de la pression de la solution optimale (CFX)

 72

4.2 Cas de l’aspirateur Porjus U9

La turbine U9 installée àla centrale hydroélectrique Porjus est une turbine Kaplan de 10 MW, sur
laquelle plus de 200 capteurs ont été installés pour la mesure des forces, de la pression, des
contraintes et des déplacements. Ces dispositifs font de cette installation un laboratoire
hydroélectrique à grande échelle, en mesure de fournir des données détaillées pour le
développement et l’amélioration des turbines hydrauliques. L’aspirateur de la turbine Porjus U9 a
été choisi comme cas test d’aspirateur coudé dans ce projet, en raison de la disponibilité des
données expérimentales.

(Mulu & Cervantes, 2009) ont décrit le banc d’essai détaillé et les résultats expérimentaux pour
l’ensemble de la turbine Porjus U9. Dans les essais sur l’aspirateur U9, les profils de vitesse ont
étémesurés sur trois sections (I, II, III) àtravers quatre fenêtres (a, b, c, d) avec la vélocimétrie
Laser Doppler (LDV).

Figure 4.19 Banc d’essai de l’aspirateur Porjus U9 (Mulu & Cervantes, 2009)

4.2.1 Maillage

DTMesh, un code de génération de maillage interne dédiéàla génération de maillages structurés

pour les aspirateurs a étéutilisépour générer les maillages pour l’aspirateur Porjus U9. Le maillage
a étégénéréen utilisant le valeurs de defaut de DTMesh. Le maillage de l’aspirateur U9 contient
683,866 nœuds et son entrée est la section I de la Figure 4.19. La longueur de la partie d’extension
située àla sortie de l’aspirateur a été fixée à 3 fois le diamètre de l’entrée.
 73

Figure 4.20 Maillage de l’aspirateur Porjus U9

4.2.2 Validation avec les données expérimentales

Dans les essais, les profils de vitesse àchaque section ont étémesurés au moyen de quatre fenêtres,
ce qui a fourni quatre profils de vitesse légèrement différents dans chaque section (Mulu &
Cervantes, 2009). Cependant, dans la simulation, la condition aux limites d’entrée des vitesses est
un profil de vitesse 1D, ce qui suppose que les profils de vitesse d’entrée sont axisymétriques. Par
conséquent, un profil de vitesse moyen des quatre profils de vitesse mesurés dans la section I a été
utilisé comme condition aux limites d’entrée de la simulation.

Dans le cas d’un l’aspirateur coudé, l’erreur devrait être observée plus clairement le long de la
direction du cintrage (fenêtres a-c) que dans la direction symétrique (fenêtres b-d). Par conséquent,
les profils de vitesse de la simulation numérique dans les sections I, II, III ont étémesurés le long
des fenêtres a-c et comparés avec leurs résultats expérimentaux correspondants.
 74

Figure 4.21 Profils de vitesse mesurés et numériques dans la section I

Figure 4.22 Profils de vitesse mesurés et numériques dans la section II

 75

Figure 4.23 Profils de vitesse mesurés et numériques dans la section III

Figure 4.24 Contour de pression (OpenFOAM)

Figure 4.25 Contour de pression (CFX)

 76

Les figures ci-dessus montrent la comparaison entre les profils de vitesse mesurés et les profils de
vitesse numériques. La Figure 4.21 montre les profils de vitesse àla section I, qui est la condition
aux limites d’entrée, et CFX semble avoir la meilleure capacité pour capturer les détails de la
condition aux limites des profils de vitesse d’entrée. Cependant, la Figure 4.23 montre que les
répartitions de vitesses prédites par OpenFOAM correspondent mieux aux données expérimentales
que celles prédites par CFX. Les Figure 4.24 et Figure 4.25 présentent les contours de pression
prédits par OpenFOAM et CFX, qui sont 30 contours variant uniformément de -20000 Pa à800 Pa.

4.2.3 Solution initiale pour l’optimisation

Figure 4.26 Courbe de facteur de perte d’énergie en fonction du nombre de tourbillon (Porjus U9)

Tel qu’illustré à la Figure 4.26, le nombre de tourbillon optimal a étéobtenu et utilisécomme point
de départ dans la méthodologie d’optimisation proposée dans ce projet (Figure 4.27, Figure 4.28 et
Figure 4.29). Au BEP, le facteur de perte d’énergie est 0,1585.
 77

Figure 4.27 Point de départ d’optimisation (Porjus U9)

Figure 4.28 Champ de vitesse du point de départ (Porjus U9)

Figure 4.29 Contour de pression du point de départ (Porjus U9)

 78

4.3 Résultats d’optimisation

4.3.1 Cas # 1

Dans le cas # 1, toutes les configurations d’optimisation sont identiques au cas du diffuseur
conique. Seul le maillage a été remplacé et les paramètres de dtpmaker ont été modifiés pour
correspondre au nouveau maillage de l’aspirateur Porjus U9. Plus précisément, l’épaisseur de la
couche aux limites extérieure dans dtpmaker a été réduite de 5% à 2% du diamètre, et l’angle
d’ouverture qui est utilisé pour le calcul du profil de vitesse radiale a étérégléà6 degrés. Compte
tenu du fait que l’écoulement turbulent est plus chaotique dans l’aspirateur Porjus U9 que dans le
diffuseur conique, la contrainte résiduelle finale de ce cas a étéaugmentée à10−3 pour conserver
un plus grand nombre de solutions potentielles.

Le graphique de convergence de ce cas d’optimisation est montréàla Figure 4.30.

Figure 4.30 Graphique de convergence de l’optimisation (Cas # 1)

La Figure 4.31 montre le profil de vitesse d’entrée optimale trouvé par NOMAD dans le cas # 1.
Les Figure 4.32 et Figure 4.33 présentent les courbes de convergence du résidu et la courbe de
facteur de perte d’énergie obtenus àpartir de la simulation de la solution optimale.
 79

Figure 4.31 Solution optimale (Cas # 1)

Figure 4.32 Courbes résiduelles de la solution optimale (Cas # 1)

 80

Figure 4.33 Courbe de facteur de perte d’énergie de la solution optimale (Cas # 1)

Le profil d’entrée optimal de vitesse a été imposé afin d’effectuer une simulation avec CFX pour
fin de comparaison. Le facteur de perte d’énergie prédit par OpenFOAM est de 0,0565451 et celui
prédit par CFX est de 0,0548490. Les solutions d’écoulement sont visualisées comme Figure 4.34,
Figure 4.35, Figure 4.36 et Figure 4.37.

Figure 4.34 Champ de vitesse de la solution optimale (Cas # 1, OpenFOAM)

 81

Figure 4.35 Champ de vitesse de la solution optimale (Cas # 1, CFX)

Figure 4.36 Contour de pression la solution optimale (Cas # 1, OpenFOAM)

Figure 4.37 Contour de pression la solution optimale (Cas # 1, CFX)

 82

4.3.2 Cas # 2

Dans le cas # 2, les contraintes d’optimisation ont été assouplies. Les limites inférieures des points
de contrôle du profil tangentiel ont étémises à-1 plutôt qu’à 0 comme dans le cas précédent (Figure
4.38). Cette modification permet àdtpmaker de générer un profil de vitesse tangentielle qui a une
direction positive sur une portion du rayon et négative sur une autre portion.

-1

Figure 4.38 Les gammes des points de contrôle tangentiels (Cas # 2)

 83

Le graphique de convergence de ce cas d’optimisation est montréàla Figure 4.39.

Figure 4.39 Graphique de convergence de l’optimisation (Cas # 2)

Figure 4.40 Solution optimale (Cas # 2)

 84

La Figure 4.40 montre le profil de vitesse d’entrée optimale trouvé par NOMAD dans le cas # 2.
Figure 4.41 et Figure 4.42 sont les courbes résiduelles et la courbe du facteur de perte d’énergie à
partir de la simulation de la solution optimale qui montrent la convergence.

Figure 4.41 Courbes résiduelles de la solution optimale (Cas # 2)

Figure 4.42 Courbe de facteur de perte d’énergie de la solution optimale (Cas # 2)

 85

Le profil d’entrée optimal de vitesse a été imposé afin d’effectuer une simulation avec CFX pour
fin de comparaison. Le facteur de perte d’énergie prédit par OpenFOAM est 0,0547851 et celui
prédit par CFX est de 0,0533725. Les solutions d’écoulement sont visualisées comme Figure 4.43,
Figure 4.44, Figure 4.45 et Figure 4.46.

Figure 4.43 Champ de vitesse de la solution optimale (Cas # 2, OpenFOAM)

Figure 4.44 Champ de vitesse de la solution optimale (Cas # 2, CFX)

Figure 4.45 Contour de pression de la solution optimale (Cas # 2, OpenFOAM)

 86

Figure 4.46 Contour de pression de la solution optimale (Cas # 2, CFX)

4.3.3 Discussion

Une comparaison des résultats produits par chacune des deux formulations du problème
d’optimisation ne permet d’observer que très peu de différences au niveau des pertes entre les deux
solutions trouvées. Dans le cas 2, une contrainte plus souple a étéimposée sur le profil de vitesse
tangentielle, ce qui devrait permettre àl’algorithme d’optimisation de trouver une solution optimale
dans un espace d’optimisation plus grand. Cependant, le facteur de perte d’énergie de la solution
optimale est de 0,0565451 dans le cas 1 et de 0,0547851 dans le cas 2, ce qui représente une
différence d’environ 3% inférieur au cas 1. La contrainte plus souple dans le cas 2 n’a donc pas
permis d’améliorer de manière significative le résultat d’optimisation. Au contraire, dans le cas 2,
le profil de vitesse tangentielle optimal est plus fluctuant, ce qui pourrait être plus difficile àréaliser
par le concepteur de turbine. Ainsi, la contrainte dans le cas 1 est considérée meilleure que le cas
2 pour ce problème d'optimisation.
 87

CHAPITRE 5 CONCLUSION ET RECOMMANDATIONS

5.1 Conclusion
Ce projet a permis d’étendre le travail de Galvan (Galván, Rubio, Pacheco, Mendoza, & Toledo,
2012; Sergio, Carlos, Pacheco, Gildardo, & Georgina, 2013) et a proposéune nouvelle approche
d’optimisation pour la résolution du problème d’optimisation des profils de vitesse d’entrée
d’aspirateurs. Cette méthodologie d’optimisation permet de déterminer les profils les plus
appropriés de vitesse d’entrée pour les aspirateurs en sortie des roues de Kaplan. L’objectif
principal de ce projet était d’examiner la combinaison du paramétrage du profil basé sur
l’interpolation et l’algorithme MADS et de mettre en œuvre cette nouvelle méthodologie avec le
code àsource ouvert OpenFOAM.

Comparé à la solution optimale du processus d’optimisation traditionnel des profils de vitesse

d’entrée de l’aspirateur, qui est aussi le point de départ de l’optimisation, cette nouvelle approche
réduit de près de 90% la perte d’énergie à l’intérieur de l’aspirateur dans le cas du diffuseur conique
et réduit de 60 % la perte d’énergie dans des cas de l’aspirateur Porjus U9. La mise en œuvre de la
méthodologie d’optimisation peut ainsi être considérée comme un succès.

Ce projet visait également àutiliser des logiciels àcodes source ouvert pour remplacer les codes
commerciaux. Deux logiciels àcodes source ouverts, àsavoir NOMAD et OpenFOAM, ont été
introduits dans ce projet pour mettre en œuvre l’algorithme d’optimisation MADS et pour résoudre
l’équation RANS avec fermeture k-epsilon. En plus de l’avantage d’éliminer les frais de licence,
un autre avantage apportépar les logiciels àcode source ouvert dans ce projet est qu'il n’y ait plus
aucune limitation dans l’utilisation des logiciels en parallèle, de telle sorte que la capacité de
grappes de calcul àhaute performance peut être pleinement exploitée et le temps de calcul est ainsi
réduit.

En termes de paramétrage du profil de vitesse d’entrée, la méthode de paramétrage basée sur

l'interpolation offre une plus grande flexibilité pour créer les profils de vitesse d’entrée. Cette
méthode de paramétrage a permis l’optimisation de l’algorithme afin de rechercher un plus large
espace d’optimisation de sorte que la meilleure solution puisse être obtenue.
 88

Dans la partie de simulation, beaucoup d’efforts ont été mis dans la calibration et la validation des
résultats obtenus avec OpenFOAM, car l’efficacité de la méthodologie d’optimisation est basée sur
la fidélitédes simulations CFD. Pour valider la précision des simulations, deux cas tests ont été
utilisés. Dans le cas du diffuseur conique, les simulations avec OpenFOAM ont étécalibrées avec
des simulations bien testées de CFX, qui est l’un des codes CFD commerciaux standard utilisés en
industrie dans le domaine des turbomachines. Dans le cas de l’aspirateur Porjus U9, les résultats
de simulation obtenus avec OpenFOAM et CFX ont étécomparés avec les données expérimentales
du projet de turbine Porjus U9. La validation a montréque la distribution de vitesse prédite par
OpenFOAM a un bon accord avec les données expérimentales, encore mieux que la prédiction de
CFX. Une conclusion peut donc être tirée que, en termes de projet d’optimisation de l’aspirateur,
OpenFOAM est capable de remplacer CFX pour sauver des frais de licence et pour améliorer
considérablement l’évolutivité de l’optimisation.

Dans le cas d’optimisation, l’algorithme MADS a aussi été introduit dans ce projet ce qui a réduit
le nombre d’évaluations de la fonction objectif. Par rapport àl’algorithme génétique multi-island
(MIGA) utilisé dans la méthodologie de Galvan, MADS a résolu un problème d’optimisation à une
échelle plus vaste, tout en utilisant un budget d’optimisation plus faible.

Par rapport à l’AG, l’algorithme MADS a une meilleure capacité àtraiter des contraintes non-
linéaires. Ainsi, deux contraintes non-linéaires ont été ajoutées. La contrainte du nombre de
tourbillon qui a été ajoutée, a éliminé les résultats les plus indésirables avant de lancer les
simulations. D’un point de vue statistique, la contrainte sur le nombre de tourbillon a éliminé
environ 80% des évaluations de la fonction objectif dans les 400 premières itérations, ce qui a
réduit de manière significative la consommation de temps global du processus d’optimisation. Par
ailleurs, la contrainte sur le résidu a éliminé les profils de vitesse d’entrée qui ne peuvent converger
dans la simulation CFD vers un état stable, ce qui n’est également pas souhaitable. Dans l’ensemble
ce projet a ainsi atteint son objectif général.

Dans ce projet, il y a encore plusieurs points qui pourraient être améliorés, mais qui n’ont pas été
terminés en raison de ressources et de temps limités.

Le logiciel cfdOpt a étéconçu pour résoudre le problème d'optimisation du profil de vitesse d'entrée
avec différents types de maillages d’aspirateurs. Au début, nous avons décidéd'utiliser 5 points de
contrôle pour manipuler le profil de vitesse axiale, 5 points de contrôle pour le profil de vitesse
 89

tangentielle. Tous les points de contrôle sont répartis uniformément le long du rayon d'entrée.
Cependant, lorsque cfdOpt résout le problème d'optimisation sur un maillage d’aspirateur sans
moyeu, l'épaisseur des couches limites intérieures est mise à 0 dans dtpmaker. Dans ce cas, le
premier point de contrôle du profil de vitesse tangentielle est situé au centre de l'entrée. Si la
variable de conception correspondante n'est pas égale àzéro, il existera deux vecteurs de vitesse
au centre de l'entrée indiquant des directions opposées, ce qui est physiquement irréaliste. Ainsi,
un traitement spécial est nécessaire pour les maillages des aspirateurs sans moyeu. Dans ce projet,
ce problème a étécontournéen fixant simplement les variables de conception correspondantes à
zéro dans NOMAD.

L’algorithme MADS a donc toujours résolu un problème d’optimisation comportant 10 variables

de conception, entrainant ainsi une augmentation du nombre d’évaluations. Ce problème n'a pas
été parfaitement résolu en raison du temps limité. L’algorithme MADS devrait avoir de meilleures
performances si cette partie avait étéajustée.

5.2 Travail futur

La solution optimale trouvée dans ce projet permet d’obtenir une faible perte d’énergie pour une
seule condition de fonctionnement de l’aspirateur. Le profil idéal de vitesse d’entrée ne devrait pas
seulement bien performer au BEP mais aussi garder la perte d’énergie à un niveau bas lorsque l’état
de fonctionnement s’en éloigne. Cependant, l’efficacité hors-BEP des solutions optimales qui ont
ététrouvées par ce projet ne sont pas encore claires. Des efforts supplémentaires devraient être mis
dans le développement d’un processus de post-optimisation qui permettrait d’évaluer l’efficacité
de la perte d’énergie de la solution optimale pour plusieurs conditions de fonctionnement. Un
procédépotentiel de post-optimisation pourrait être similaire au test de rotation de corps solide, le
profil de vitesse axiale est fixe et le profil de vitesse tangentielle est ajustépour obtenir un nombre
de tourbillon différent.

Avec la méthode de paramétrage basésur l’interpolation il existait aussi un potentiel d’amélioration.

Dans la configuration actuelle, le profil de vitesse axiale et le profil de vitesse tangentielle ont le
même nombre de points de contrôle et tous les points de contrôle ont été simplement répartis
uniformément le long du rayon. Cependant, le résultat a montréque la variation du profil de vitesse
n’est pas uniforme. Il y une plus grande variation dans la partie extérieure que dans la partie
 90

intérieure, et le profil de vitesse tangentielle présente une plus grande variation le long du rayon
que le profil de vitesse axiale. Par conséquent, plus de points de contrôle devraient être alloués à
ces zones pour un réglage plus fin.
 91

BIBLIOGRAPHIE

Abramson, M. A., Audet, C., Dennis Jr, J. E., & Digabel, S. L. (2009). OrthoMADS: A
deterministic MADS instance with orthogonal directions. SIAM Journal on optimization,
20(2), 948-966.

Andersson, U., & Dahlbäck, N. (1998). Experimental evaluation of draft tube flow—a test case for
CFD simulations. Paper presented at the XIX IAHR Symposium on Hydraulic Machinery
and Systems, Singapore.

Ansys, C. (2012). ANSYS CFX-solver theory guide. ANSYS CFX Release, 11, 69-118.

Argyropoulos, C., & Markatos, N. (2015). Recent advances on the numerical modelling of
turbulent flows. Applied Mathematical Modelling, 39(2), 693-732.

Audet, C., & Dennis Jr, J. E. (2002). Analysis of generalized pattern searches. SIAM Journal on
optimization, 13(3), 889-903.

Audet, C., & Dennis Jr, J. E. (2006). Mesh adaptive direct search algorithms for constrained
optimization. SIAM Journal on optimization, 17(1), 188-217.

Bahrami, S., Tribes, C., Devals, C., Vu, T. C., & Guibault, F. (2013). Multi-Objective Optimization
of Runner Blades Using a Multi-Fidelity Algorithm. Paper presented at the ASME 2013
Power Conference.

Bahrami, S., Tribes, C., von Fellenberg, S., Vu, T., & Guibault, F. (2014). Multi-fidelity design
optimization of Francis turbine runner blades. Paper presented at the IOP Conference
Series: Earth and Environmental Science.

Bahrami, S., Tribes, C., von Fellenberg, S., Vu, T. C., & Guibault, F. (2015). Physics-based
Surrogate Optimization of Francis Turbine Runner Blades, Using Mesh Adaptive Direct
Search and Evolutionary Algorithms. International journal of fluid machinery and systems,
8(3), 209-219.

Chow, E., Falgout, R. D., Hu, J. J., Tuminaro, R. S., & Yang, U. M. (2006). A survey of
parallelization techniques for multigrid solvers. Parallel processing for scientific
computing, 20, 179-201.

Clarke, F. H. (1990). Optimization and nonsmooth analysis (Vol. 5): Siam.

 92

Conn, A. R., Scheinberg, K., & Vicente, L. N. (2009). Introduction to derivative-free optimization
(Vol. 8): Siam.

Fritsch, F. N., & Carlson, R. E. (1980). Monotone piecewise cubic interpolation. SIAM Journal on
Numerical Analysis, 17(2), 238-246.

Gagnon, J.-M., Aeschlimann, V., Houde, S., Flemming, F., Coulson, S., & Deschenes, C. (2012).
Experimental investigation of draft tube inlet velocity field of a propeller turbine. Journal
of Fluids Engineering, 134(10), 101102.

Gagnon, J., Iliescu, M., Ciocan, G., & Deschênes, C. (2008). Experimental investigation of runner
outlet flow in axial turbine with LDV and stereoscopic PIV. Paper presented at the 24th
IAHR Symposium on Hydraulic Machinery and Systems, Foz do Iguassu, Brazil, October.

Galván, S., Reggio, M., & Guibault, F. (2011). Assessment Study of K-ɛ Turbulence Models and
Near-Wall Modeling for Steady State Swirling Flow Analysis in Draft Tube Using Fluent.
Engineering Applications of Computational Fluid Mechanics, 5(4), 459-478.

Galván, S., Reggio, M., & Guibault, F. (2012). Optimization of the inlet velocity profile in a conical
diffuser. Paper presented at the ASME 2012 Fluids Engineering Division Summer Meeting
collocated with the ASME 2012 Heat Transfer Summer Conference and the ASME 2012
10th International Conference on Nanochannels, Microchannels, and Minichannels.

Galván, S., Rubio, C., Pacheco, J., Mendoza, C., & Toledo, M. (2012). Optimization methodology
assessment for the inlet velocity profile of a hydraulic turbine draft tube: part I—computer
optimization techniques. Journal of Global Optimization, 55(1), 53-72.
doi:10.1007/s10898-012-9946-8

Galvan, S. R. (2007). Optimization of the inlet velocity profile of the Turbine 99 draft tube.
(NR35512 Ph.D.), Ecole Polytechnique, Montreal (Canada), Ann Arbor. Retrieved from
http://search.proquest.com/docview/304718911?accountid=12339

http://mcgill.on.worldcat.org/atoztitles/link?sid=ProQ:&issn=&volume=&issue=&title=Optimiza
tion+of+the+inlet+velocity+profile+of+the+Turbine+99+draft+tube&spage=&date=2007
-01-
01&atitle=Optimization+of+the+inlet+velocity+profile+of+the+Turbine+99+draft+tube&
 93

au=Galvan%2C+Sergio+Ricardo&id=doi: ProQuest Dissertations & Theses Full Text

database.

Golberg, D. E. (1989). Genetic algorithms in search, optimization, and machine learning. Addion
wesley, 1989, 102.

Grotjans, H. (2001). Simulations on draft tube with CFX. Paper presented at the Proceedings of
Turbine-99, Second ERCOFTAC Workshop on Draft Tube Flow, Vattenfall Utveckling
AB, Älvkarleby, Sweden. .

Keck, H., & Sick, M. (2008). Thirty years of numerical flow simulation in hydraulic
turbomachines. Acta Mechanica, 201(1-4), 211-229.

Launder, B. E., & Spalding, D. (1974). The numerical computation of turbulent flows. Computer
methods in applied mechanics and engineering, 3(2), 269-289.

Liu, S., Wu, Y., Wu, X., & Nishi, M. (2011). Optimization of a Francis Turbine Through Flow
Simulation with a RNG k-ω Turbulence Model. International Journal of Turbo and Jet
Engines, 28(1), 53-58.

Marjavaara, B., & Lundström, T. (2006). Redesign of a sharp heel draft tube by a validated CFD‐
optimization. International Journal for Numerical Methods in Fluids, 50(8), 911-924.

Moler, C. B. (2008). Numerical Computing with MATLAB: Revised Reprint: Siam.

Mulu, B., & Cervantes, M. (2009). Experimental investigation of a Kaplan model with LDA. Paper
presented at the Proceedings of the Water Engineering for a Sustainable Environment
Congress: 33rd IAHR Congress, Vancouver, Canada.

Nilsson, H., & Cervantes, M. J. (2012). Effects of inlet boundary conditions, on the computed flow
in the Turbine-99 draft tube, using OpenFOAM and CFX. IOP Conference Series: Earth
and Environmental Science, 15(3), 032002.

Nocedal, J., & Wright, S. (2006). Numerical optimization: Springer Science & Business Media.

P Skåre, O. D. (2001). CFD calculations of the flow in Kaplan draft tube. Paper presented at the
Proceedings of Turbine-99, Second ERCOFTAC Workshop on Draft Tube Flow,
Vattenfall Utveckling AB, Älvkarleby, Sweden.
 94

Petit, O., Nilsson, H., Vu, T., Manole, O., & Leonsson, S. (2008). The flow in the U9 Kaplan
turbine—preliminary and planned simulations using CFX and OpenFOAM. Paper
presented at the Proceedings of the 24th IAHR Symposium on Hydraulic Machinery and
Systems.

Puente, L., Reggio, M., & Guibault, F. (2003). Automatic shape optimization of a hydraulic turbine
draft tube. Paper presented at the Proceedings of the International Conference CFD2003,
Vancouver, BC.

Sergio, G., Carlos, R., Pacheco, J., Gildardo, S., & Georgina, C. (2013). Optimization methodology
assessment for the inlet velocity profile of a hydraulic turbine draft tube: part II—
performance evaluation of draft tube model. Journal of Global Optimization, 55(4), 729-
749.

Shih, T.-H., Liou, W., Shabbir, A., Yang, Z., & Zhu, J. (1994). A new k-epsilon eddy viscosity
model for high Reynolds number turbulent flows: Model development and validation.

Susan-Resiga, R., Ciocan, G. D., Anton, I., & Avellan, F. (2006). Analysis of the swirling flow
downstream a Francis turbine runner. Journal of Fluids Engineering, 128(1), 177-189.

Susan-Resiga, R., Muntean, S., Avellan, F., & Anton, I. (2011). Mathematical modelling of
swirling flow in hydraulic turbines for the full operating range. Applied Mathematical
Modelling, 35(10), 4759-4773.

Susan-Resiga, R., Muntean, S., Ciocan, T., Joubarne, E., Leroy, P., & Bornard, L. (2012). Influence
of the velocity field at the inlet of a Francis turbine draft tube on performance over an
operating range. Paper presented at the IOP Conference Series: Earth and Environmental
Science.

Van Dyke, B., & Asaki, T. J. (2013). Using QR decomposition to obtain a new instance of mesh
adaptive direct search with uniformly distributed polling directions. Journal of
Optimization Theory and Applications, 159(3), 805-821.

Vu, T. C., Devals, C., Zhang, Y., Nennemann, B., & Guibault, F. (2011). Steady and unsteady flow
computation in an elbow draft tube with experimental validation. International journal of
fluid machinery and systems, 4(1), 85-96.
 95

Wetter, M., & Wright, J. (2003). Comparison of a generalized pattern search and a genetic
algorithm optimization method. Paper presented at the Proceedings of the 8th International
IBPSA Conference, Eindhoven, Netherlands.

Wilcox, D. C. (1988). Reassessment of the scale-determining equation for advanced turbulence

models. AIAA journal, 26(11), 1299-1310.

Wilcox, D. C. (2008). Formulation of the kw turbulence model revisited. AIAA journal, 46(11),
2823-2838.
 96

ANNEXE A – L’EXPRESSION DU POST-TRAITEMENT

OpenFOAM (diffuseur conique, surface)
FoamFile {
version 2.0;
format ascii;
class dictionary;
object controlDict;
}
functions {
outlet_surface {
type createSampledSurface;
surfaceName diffuserOutlet;
surface {
type plane;
basePoint (0 0 2.3523);
normalVector (0 0 1);
interpolate true;
}
outputControl timeStep;
outputInterval 1;
}
energyLossFactor {
type swakExpression;
variables (
"Vout{surface'diffuserOutlet} = sum(Sf()&U)/(sum(area()));"
"Pout{surface'diffuserOutlet} = sum(p * area())/sum(area());"
"Vthroat{inlet} = sum(phi) / (sum(area()));"
"Athroat{inlet} = sum(area());"
);
valueType patch;
patchName inlet;
outputControlMode timeStep;
outputInterval 1;
accumulations (
sum
);
expression "((p+0.5*magSqr(U)) - (Pout + 0.5 * Vout * Vout)) * area() / (0.5 * Vthroat * Vthroat *
Athroat)";
verbose true;
}
}
 97

OpenFOAM (diffuseur conique, faceZone)

FoamFile {
version 2.0;
format ascii;
class dictionary;
object faceSetDict;
}
name diffuserOutlet;
action new;
topoSetSources (
expressionToFace {
expression "fpos().z > interpolate(2.3522) && fpos().z < interpolate(2.3524)";
}
);
FoamFile {
version 2.0;
format ascii;
class dictionary;
object controlDict;
}
functions {
energyLossFactor {
type swakExpression;
variables (
"Vout{faceZone'diffuserOutlet} = sum(phi * flip())/(sum(area()));"
"Pout{faceZone'diffuserOutlet} = sum(p * area())/sum(area());"
"Vthroat{inlet} = sum(phi) / (sum(area()));"
"Athroat{inlet} = sum(area());"
);
valueType patch;
patchName inlet;
outputControlMode timeStep;
outputInterval 1;
accumulations (
sum
);
expression "((p+0.5*magSqr(U)) - (Pout + 0.5 * Vout * Vout)) * area() / (0.5 * Vthroat * Vthroat *
Athroat)";
verbose true;
}
}

CFX (diffuseur conique)

 98

LIBRARY:
CEL:
EXPRESSIONS:
AreaOut = area()@diffuserOutlet
Athroat = area()@Inlet
DensityExpr = 997 [kg m^-3]
DthroatExpr = maxInletRadius * 2
FluxIn = areaInt(Velocity w)@Inlet
FluxOut = areaInt(Velocity w)@diffuserOutlet
Kin = areaAve(0.5*DensityExpr*Velocity^2)@Inlet
KineticEnergyExpr = 0.5*DensityExpr*VelocityExpr*VelocityExpr
Kout = maxVal(0.5 * DensityExpr * Vout*Vout)@diffuserOutlet
Kthroat = maxVal(0.5 * DensityExpr * Vthroat*Vthroat)@Inlet
LossDT = (Ptin - PtoutLab)/Kthroat
Pstatic = PstaticExpr
PstaticExpr = (Pressure)*1
Ptin = areaAve(Pressure + 0.5 * DensityExpr * Velocity * Velocity )@Inlet
PtotalExpr = Pstatic + KineticEnergyExpr
PtoutLab = areaAve(PstaticExpr)@diffuserOutlet + Kout
VelocityExpr = sqrt(Velocity u*Velocity u + Velocity v*Velocity v + Velocity w*Velocity w)
Viscosity ratio = Eddy Viscosity/Dynamic Viscosity
Vout = mFluxOut/(AreaOut*DensityExpr )
Vthroat = mFluxIn/(Athroat*DensityExpr)
inletP = areaAve(Pressure/DensityExpr)@Inlet
mFluxIn = areaInt(DensityExpr * Velocity w)@Inlet
mFluxOut = areaInt(DensityExpr * Velocity w)@diffuserOutlet
maxInletRadius = maxVal(sqrt(x^2 + y^2))@inletRadius
outletP = areaAve(Pressure/DensityExpr)@diffuserOutlet
swirl = lengthInt(Density * Velocity w * (x^2 + y^2) * sqrt(Velocity u^2 + Velocity
v^2) )@inletRadius / (maxInletRadius * lengthInt(Density * sqrt(x^2 + y^2) * Velocity
w^2 )@inletRadius )
END
END
END